I am trying to design a linear-phase FIR filter with the following frequency and phase responses:
Designing the filter
Given its characteristics (peak filter at 1kHz with 9dB gain and Q=7), I got the following IIR biquad coefficients:
$$b_0 = 1.0889938082113808$$ $$b_1 = -1.8806408829603927$$ $$b_2 = 0.80787505079350063$$ $$a_0 = 1$$ $$a_1 = -1.8806408829603927$$ $$a_2 = 0.89686885900488144]$$
I then computed the impulse response of that biquad filter and truncated it to 1024 taps. I got the following filter:
To make the phase of my FIR filter linear, I applied an FFT on its impulse response, I put the imaginary part to 0 and came back to time domain by doing an inverse FFT.
The final taps I got can be found there.
Now, when I plot the frequency/phase responses of this filter using Numpy/Scipy or Matlab, I get exactly the same plot as the one of my expected filter (see figure 1). But, when I apply this filter on a pink noise signal by convolution and measure the transfer function, I get this instead:
I am pretty lost here, I would assume that convolving my noise signal with my filter and measuring the transfer function would give me the same curve as the one I computed theoritically but it seems not to be the case (that is, in my mind, the definition of what convolution does).
Some additional information:
- I have verified that the problem does not come from my implementation of the convolution operation, by testing with many different implementations; if you still believe the problem come from there, you can try yourself with the taps I provided.
- I trust the software I used to measure the transfer function depicted in the last plot, and anyway I could hear that the transfer function is not the one I expect.
To plot the transfer functions from Python, I used the following code:
import matplotlib.pyplot as plt import numpy as np from scipy import signal def plot_transfer_function(b, a=1): """Plot the frequency and phase response of a digital filter.""" w, h = signal.freqz(b, a) fig = plt.figure() plt.title('Digital filter frequency response') # Frequency response ax1 = fig.add_subplot(211) plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), 20 * np.log10(abs(h)), 'b') # Make the plot pretty plt.ylabel('Amplitude [dB]', color='b') plt.xlabel('Frequency [Hz]') plt.axis([20, 20000, -12, 12]) plt.grid() # show the grid # Phase response ax2 = fig.add_subplot(212) angles = np.unwrap(np.angle(h)) plt.semilogx((w * SAMPLING_RATE) / (2 * np.pi), np.degrees(angles), '#aa0000') # Make the plot pretty plt.ylabel('Angle (degrees)', color='#aa0000') plt.axis([20, 20000, -180, 180]) plt.grid() fig.axes.set_axis_off() # disable main figure axis plt.show() plot_transfer_function(TAPS) # "TAPS" is an array containing the 1024 taps
So, what is wrong with my reasoning here? Why do I not measure the same frequency/phase response when calculating it and when measuring it after a convolution?
Additional question: for my FIR filter to have linear phase, I thought it should be symmetric (or anti-symmetric). Well, when I plot my 1024 taps, they appear to be pretty symmetric, excepted the first which does not have its equivalent at the end of the impulse response. Why is it so?