1
$\begingroup$

I am trying to calculate an actual frequency response from a desired frequency response.

fftSize = 128;
fs = 8000;

I have 2 input arrays: f[] and aDesired[], both of which are Real and length fftSize+1 (DC + Nyquist included). The f[] array is linearly spaced between 0 and fs/2 (0:4000). The aDesired[] array is either gains in dB or linear values.

I can easily change the length of f[] and aDesired[] as I have an interpolation function which gives any size outputs, but i thought this was a good number.

I would like to take an IFFT of ADesired[] to generate coefficients (with a configurable length), window and truncate the coeffs, and then FFT to show what the actual response would be, aActual[].

How can I take an IFFT/FFT of this data when the inputs are both Real? Or how do I prepare the data more?

Finally, it doesnt matter about the phase, if its linear or not it doesnt matter.

Is there maybe a library that could be of use with my input data?

Summary:

Input: f[129] : linearly spaced fc's
       aDesired[129] : Real values of linearly spaced gains at fc's (either in lin or log),

Processing:
       coeffs = IFFT(f, aDesired)
       coeff_new Truncate(window * coeff)
       aActual = FFT(f, coeff_new)

Output: aActual[]

Many thanks

$\endgroup$
  • 3
    $\begingroup$ What you are describing is essentially the frequency sampling method for designing arbitrary FIR filters. $\endgroup$ – Jason R Jun 11 '13 at 2:24
3
$\begingroup$

As Jason R pointed out in his comment, your way of obtaining an actual frequency response is just one way of doing it. There are many other ways of determining an actual filter response. Any other frequency domain FIR filter design method could do the job, and with each method you would get a slightly different result. Having said that, let's have a look at the way you want to do it:

You have a real-valued desired response. This is the magnitude of the desired frequency response, because for causal (realizable) systems, the frequency response is always complex-valued. However, this is no problem since any finite length impulse response can always be made causal by adding some delay. First, we need to make the desired response symmetric, as is necessary for real-valued systems. Let the column vector D be the vector of desired magnitude values on the equidistant frequency grid between 0 and $f_s/2$, then the symmetric desired response is (using Matlab/Octave syntax):

Dsym=[D;D(N/2:-1:2)];

where $N$ is the FFT length. We can obtain the (non-causal) impulse response by doing an inverse FFT:

h=real(ifft(Dsym));

Note that taking the real part is only necessary to remove the tiny imaginary part due to numerical errors. To make the impulse response causal, we need to shift the second half of it to the beginning:

h_c=[h(N/2+2:N); h(1:N/2+1)];

The impulse response can now be multiplied with a window function:

h_w = h_c.*w;

Your actual frequency response is then

H=fft(h_w);

It is complex valued, so in order to compare it to the desired response, you need to plot its magnitude:

plot(abs(H))

The phase of the actual response is linear, i.e. the system just adds some delay but does not introduce any phase distortions.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hi Matt Thanks very much for your answer, its exactly what i need, however i haven't been able to implement it as I would expect. I think maybe my numbers are not quite correct, because the actual and desired are quite different. My input, D, is 129 samples with DC and Nyquist, is this what you expect? $\endgroup$ – Andy Jun 11 '13 at 19:03
  • $\begingroup$ It depends a lot on your choice of D. If you add the definition of D to your question then I can take a look. $\endgroup$ – Matt L. Jun 11 '13 at 19:06
  • $\begingroup$ N = 128; D = zeros(129,1); D(20:100) = ones*3; D(50:75) = ones*6; Dsym=[D;D(N:-1:2)]; % Changed from Dsym=[D;D(N/2:-1:2)] which gave length 192, is that right? h=real(ifft(Dsym,N)); h_c=[h(N/2+2:N); h(1:N/2+1)]; w = hamming(N); h_w = h_c.*w; H=fft(h_w); $\endgroup$ – Andy Jun 11 '13 at 19:16
  • $\begingroup$ sorry the formatting is not so nice, here is how i do the plots: figure; subplot(2,1,1), plot(D), TITLE('D'), axis([0 256 0 7]) subplot(2,1,2), plot(Dsym), TITLE('Dsym=[D;D(N:-1:2)] 256x1: missing final point??'), axis([0 256 0 7]) figure; subplot(2,2,1), plot(h), TITLE('h'), axis([0 128 -2 4]) subplot(2,2,2), plot(h_c), TITLE('h_c'), axis([0 128 -2 4]) subplot(2,2,3), plot(w), TITLE('w'), axis([0 128 0 1]) subplot(2,2,4), plot(h_w), TITLE('h_'), axis([0 128 -2 4]) figure; plot(D, 'r'), axis([0 130 0 7]) hold plot(abs(H), 'g'), TITLE('D vs A' ), axis([0 130 0 7]) $\endgroup$ – Andy Jun 11 '13 at 19:25
  • $\begingroup$ If you have 129 frequency points between 0 and Nyquist, then your FFT length must be N=256, i.e. then the N/2 for Dsym should be correct. The FFT also contains the redundant points between fs/2 and fs (or negative frequencies, if you like). $\endgroup$ – Matt L. Jun 11 '13 at 19:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.