Butterworth filtering behaving unexpectedly -MATLAB

Question

I'm new to signal processing but I'm trying to understand how the Butterworth filter works. To do so, I did the following: I assumed a sampling frequency of 100Hz and wanted to attentuate everything after half the Nyquist frequency. The code I used is the following

[b,a] = butter(2,0.5);
fs = 100; %Sampling frequency
t = 0:1/fs:1-1/fs;
y=0.123*sin(2*pi*20*t+0.234)+0.123*sin(2*pi*10*t); %Pair of sine curves
plot(abs(fft(y))) %This behaves as expected, with two peaks at 10 and 20
dataOut = filter(b, a, y);
figure,
plot(abs(fft(dataOut))) 

I expected the second plot to be identical to the first because there are no frequencies above the Butterworth filter's cutoff anyway. However, the signal is artifically inflated (by a small amount) for all components below 25Hz and then goes to zero. See the original signal and the filtered signal below.

Is this normal behavior? If so, it is effectively adding noise in my passband, is it not? Also, this effect is a lot worse if I use a higher order Butterworth filter.

Answer 1

The effect you observed is mainly due to spectral leakage. It is a property of the DFT and it is only indirectly related to the Butterworth filter. Note that there are an integer number of periods of both sinusoidal components of y inside the DFT window. That's why you see two nice peaks in the frequency domain without any leakage. If you truncate y only by a few samples, you'd get spectral leakage even without the Butterworth filter, because now the period of the signal doesn't line up anymore with the DFT window:

plot(abs(fft(y(1:97))))

The output signal of the Butterworth filter shows some transient before it converges to its steady state. So it is not purely sinusoidal, and - more importantly - the number of periods of the strong sinusoidal components are not aligned with the length of the DFT window. If you remove the initial transient and if you make sure the DFT window aligns with the periods of the sinusoids you get two perfect spectral peaks:

plot(abs(fft(dataOut(11:100))))

Answer 2

Butterworth filters have a flat magnitude response in the passband, which is often referred to as a maximally-flat magnitude response. The passband of a lowpass filter is typically the frequency where the magnitude response is 3-dB "down" or -3 dB of the magnitude in the passband. This frequency is often referred to as the cutoff frequency of the lowpass filter.

You can examine the frequency response of your filter using the syntax "freqz(b,a,1024,100)", where b and a are the Butterworth filter coefficients, 1024 is the number of points where the frequency response is computed (between 0 Hz and the Nyquist frequency; here the Nyquist frequency is 50 Hz). The number of points (i.e. 1024 in this case) is arbitrary.

You've designed a second-order filter. As a result, the attenuation of the magnitude response will be more gradual than for a filter with a higher order. This means frequencies below the cutoff frequency for a lower order filter will be attenuated more than for a higher order filter. See the two frequency responses below, one for the second-order filter and the other for an eighth-order filter.

If you want to know the specific amount of attenuation that the filter applies to each of your tone frequencies (i.e. 20 Hz and 10 Hz, respectively), you could filter each tone component of your composite signal and determine it's amplitude in the steady-state portion of the filter's response. Keep in mind that you've define the amplitude to be 0.123 for both tones.

One final note for the sake of clarification. The normalized cutoff you've used is 0.5 pi radians/sample. This means the filter's cutoff frequency in Hertz will be 0.5 * 50 Hz -> 25 Hz.

The FFT is a great numerical tool, but it also has many pitfalls. Scaling and spectral leakage are likely the two most common ones.