Minimum phase transformation

[takoha]

Hello all,
I am new in this forum, I begin in audio processing and to top it all, I am French: sorry for my english.

I try to convert FIR filter with linear phase, into one with the same magnitude but minimum phase.
I made an algo from this document, in the hilbert transform description
http://www.hpl.hp.com/personal/Niranjan ... nPhase.pdf

my algo is from this formula (may be misinterpreted):
teta = -j* DFT ( sign * IDFT (log (DFT(x[n])) ) (=atan (Im/Re) of the complex array result)
Then i take the complexe valu formed by the initial magnitude and this angle in the polar form, and make a IFFT to get impulse response.

First try with a simple Windowed-Sinc highpass FIR, the phase reponse look good, but seems to be scaled less than 25 times the value of a comparable IIR filter (approximatly).
After applying this factor, i get a perfect minimum phase version of the linear phase filter: exact magnitude and all left coef =0.
Second try, with a dirac with a random linear phase equalization.
Now i get a impulse with the same amplitude but not completly minimum phase: left coef are not all zero.
Phase seems to have been delayed but not enough.
It also appears that it does not depend on the previous coefficient.

Here some picture that tell best than my english
Good
http://avielbaz.free.fr/images_forum/fir1.png
http://avielbaz.free.fr/images_forum/ht1.png

Problem
http://avielbaz.free.fr/images_forum/fir2.png
http://avielbaz.free.fr/images_forum/ht2.png
(sorry for the url writting it doesn t work before 5 post)

Since i am realy a beginner, this problem may make you laugh.
but...

Cheers,
Avi

[takoha]

up

[Aleksey Vaneev]

Here's how I'm doing mp transform (this does not require H(z)^2 as some methods):

: N = kernelsize (power of 2)
: fft.size = N
: Kernel: buffer of size N
: tmp: buffer of size N / 2 + 1 - used to store original spectrum magnitude
: fftlog, fftexp - functions taking complex log() and exp() in-place
: mphom - homomorphic function (multiplies second half of the cepstrum buffer by -1)

***
fft.forward( Kernel );
fftlog( Kernel, tmp, N );
fft.inverse( Kernel );

mphom( Kernel, N );

fft.forward( Kernel );
fftexp( Kernel, tmp, N );
fft.inverse( Kernel );
***

you'll obviously need to normalize the output depending on your FFT routine which may introduce 1/N gain, for example.

Also, one half of the Kernel should be filled with zeroes or otherwise some aliasing may occur.

[takoha]

Hello Aleksey,
many thanks for the algo
I tried it, there must be things that I misunderstood (maybe because i use complex fft)
here is the code i implemented:

//kernel is 2^n samples half filled with zero of a linear phase FIR filter

//fft.forward( Kernel ) :
Complex[] temp=fft.forward( Kernel ); // store result complex in temp, with reel part in[0,N/2], symetric conjugate in [N/2+1, N-1]
//fftlog( Kernel, tmp, N ):
for(i=0;i<N;i++)
temp= log(temp); //log is the complex log function
//fft.inverse( Kernel ):
Kernel=fft.inverse( temp );

//mphom( Kernel, N );
for(i=N/2;i<N;i++)
Kernel*= -1;

//fft.forward( Kernel ):
temp=fft.forward( Kernel );
//fftexp( Kernel, tmp, N ):
for(i=0;i<N;i++)
temp= exp(temp); //exp is the complex exp function
//fft.inverse( Kernel ):
Kernel=fft.inverse( temp );

result here http://avielbaz.free.fr/images_forum/ht3.png
i checked and checked again, i don't see the problem
Or maybe it s not the good way to solve my problem (convert a linear phase FIR filter into a IIR like one with minimum group delay, i want to use them with a vst convolver plugin) and so there is something about the theory that I did not understand.

Cheers,
Avi

[Aleksey Vaneev]

First of all, I use real filter and FFT for real numbers. Beside that complex data after FFT is stored in pairs.

You've also missed the point about storing the actual magnitude of the spectrum. I'm doing the following: taking a=sqrt(r^2+i^2), storing it into the temp variable, and then taking log(a), and storing it to the Kernel, with imaginary part zeroed. So, in the essence you should not write Kernel=fft.inverse( temp ); but should write Kernel=fft.inverse( Kernel ); so that magnitudes are preserved: they are used when exponent is taken.

[takoha]

Hello Aleksey,

Oki, so I have a=log(sqrt(r^2+i^2))
I make a inverse fft of the complex value formed by re=a and im=0
Multiply half by -1
and the point i am not sur:
fftexp( Kernel, tmp, N);
Here kernel is formed by complex value, so i take the complex exp?

Avi

[Aleksey Vaneev]

and the point i am not sur:
fftexp( Kernel, tmp, N);
Here kernel is formed by complex value, so i take the complex exp?

You'll need to throw real parts from Kernel away, and use imaginary part as phase:

Kernel[i * 2] := cos (Kernel[i * 2 + 1]) * tmp ;
Kernel[i * 2 + 1] := sin (Kernel[i * 2 + 1]) * tmp ;

[takoha]

Hi Aleksey,

A little late (i stopped the Development for a while), thank you for your help
It works fine with this algo (which is close to yours)
output = ifft ( exp ( fft ( wn* ( ifft ( log|fft(input)| )))));
with wn[0]=wn[N/2]=1, wn[n=1:n=N/2-1]=2, wn[n=N/2:N-1]=0
exp is the complex exp and log real

Regards,
Avi

[Aleksey Vaneev]

Hi Aleksey,

A little late (i stopped the Development for a while), thank you for your help
It works fine with this algo (which is close to yours)
output = ifft ( exp ( fft ( wn* ( ifft ( log|fft(input)| )))));
with wn[0]=wn[N/2]=1, wn[n=1:n=N/2-1]=2, wn[n=N/2:N-1]=0
exp is the complex exp and log real

Regards,
Avi

Yep, it's the same thing except that algo I'm using retains amplitude response between the calls. Hope your filter's amplitude response does not come out squared as well.

[takoha]

Yes i do not retains the amplitude response, it's not a pb since i just want the impulse filter

Something different, is there a way to make config file for your spristine space plugin?
I want to generate it automatically (as an option) in order to make life simpler for the user (my soft is in java)

Avi

[mystran]

Can I make a really stupid question?

Why minimum phase FIR?

[Zaphod (giancarlo)]

we did it using the hilbert transformate.

[Aleksey Vaneev]

Something different, is there a way to make config file for your spristine space plugin?

It's not really possible. But you can get around that by defining presets with generic filenames like 'impulse1.wav', 'impulse2.wav', etc and then replacing these WAV files. Not the best way, of course.

[takoha]

Hi mystran

Why minimum phase FIR?

Do you mean, why minimum phase FIR vs IIR (considering that minimum phase filter can also be implemented with reccursive equation) or why minimum phase vs linear phase filter (as a matter of latency and psychoacoustics)?

Avi

[Aleksey Vaneev]

we did it using the hilbert transformate.

Could you show the algo for that?