I am trying to create an application for calculating coefficients for a graphic equalizer FIR filter. I am doing some prototyping in Matlab but I have some problems.
I have started with the following Matlab code:
% binamps vector holds 2^13 = 8192 bins of desired amplitude values for frequencies in range 0.001 .. 22050 Hz (half of samplerate 44100 Hz)
% it looks just fine, when I use Matlab plot() function
% now I get ifft
n = size(binamps,1);
iff = ifft(binamps, n);
coeffs = real(iff); % throw away the imaginary part, because FIR module will not use it anyway
But when I do the fft() of the coefficients, I see that the frequencies are stretched 2 times and the ending of my AFR data is lost:
p = fft(coeffs, n); % take the fourier transform of coefficients for a test
nUniquePts = ceil((n+1)/2);
p = p(1:nUniquePts); % select just the first half since the second half
% is a mirror image of the first
p = abs(p); % take the absolute value, or the magnitude
p = p/n; % scale by the number of points so that
% the magnitude does not depend on the length
% of the signal or on its sampling frequency
p = p.^2; % square it to get the power
sampFreq = 44100;
freqArray = (0:nUniquePts-1) * (sampFreq / n); % create the frequency array
semilogx(freqArray, 10*log10(p))
axis([10, 30000 -Inf Inf])
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
So I guess, I am using ifft wrong. Do I need to make my binamps vector twice as long and create a mirror in the second part of it? If it is the case, then is it just a Matlab's implementation of ifft or also other C/C++ FFT libraries (especially Ooura FFT) need mirrored data for inverse FFT?
Is there anything else I should know to get the FIR coefficients out of ifft?
Your frequency domain vector needs to be complex rather than just real, and it needs to be symmetric about the mid point in order to get a purely real time domain signal. Set the real parts to your desired magnitude values and set the imaginary parts to zero. The real parts need to have even symmetry such that A[N - i] = A[i]
(A[0]
and A[N / 2]
are "special", being the DC and Nyquist components - just set these to zero.)
The above applies to any general purpose complex-to-complex FFT/IFFT, not just MATLAB's implementation.
Note that if you're trying to design a time domain filter with an arbitrary frequency response then you'll need to do some windowing in the frequency domain first. You might find this article helpful - it talks about arbitrary FIR filter design usign MATLAB, in particular fir2.