Minimum phase transformation
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
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
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
- KVRAF
- 4021 posts since 7 Sep, 2002
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.
: 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.
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
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
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
- KVRAF
- 4021 posts since 7 Sep, 2002
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.
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.
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
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
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
- KVRAF
- 4021 posts since 7 Sep, 2002
You'll need to throw real parts from Kernel away, and use imaginary part as phase:takoha wrote: and the point i am not sur:
fftexp( Kernel, tmp, N);
Here kernel is formed by complex value, so i take the complex exp?
Kernel[i * 2] := cos (Kernel[i * 2 + 1]) * tmp ;
Kernel[i * 2 + 1] := sin (Kernel[i * 2 + 1]) * tmp ;
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
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
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
- KVRAF
- 4021 posts since 7 Sep, 2002
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 wrote: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
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
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
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
- KVRAF
- 7893 posts since 12 Feb, 2006 from Helsinki, Finland
-
Zaphod (giancarlo) Zaphod (giancarlo) https://www.kvraudio.com/forum/memberlist.php?mode=viewprofile&u=111268
- KVRAF
- 2596 posts since 23 Jun, 2006
we did it using the hilbert transformate.
- KVRAF
- 4021 posts since 7 Sep, 2002
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 wrote:Something different, is there a way to make config file for your spristine space plugin?
-
- KVRer
- Topic Starter
- 12 posts since 15 Nov, 2007
Hi mystran
Avi
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)?Why minimum phase FIR?
Avi
- KVRAF
- 4021 posts since 7 Sep, 2002