Matlab: convolution with bandpass filter does not cut the unwanted frequencies

Question

I have a 6ms-long signal with three frequency components sampled at 60kHz:

fs = 60000;
T = 0.006;
t = 0:1/fs:T;
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t); 

I have a bandpass filter with impulse response being the difference between two sinc functions:

M = 151;
N = 303;
n = 0:(N-1);
h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);
h(n==M) = 0.2094;

I designed a function that convolves the input with filter:

function y = fir_filter(h,x)
y = zeros(1,length(x)+length(h)-1);
for i = 1:length(x)
for j = 1:length(h)
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);
end
end

And then applied the filter:

y = fir_filter(h,x);

This produced strange results:

figure(21)
ax1 = subplot(311);
plot(x);
title('Input Signal');
ax2 = subplot(312);
plot(h);
title('FIR');
ax3 = subplot(313);
plot(y);
title('Output Signal');
linkaxes([ax1,ax2,ax3],'x')
ax2.XLim = [0,length(y)];

Since the filter is bandpass, only one frequency component was expected to survive.

I tried to use yy = filter(h,1,[x,zeros(1,length(h)-1)]); and yyy = conv(h,x); and got the same results. Please, can anybody explain me what I am doing wrong? Thank you!

Answer 1

Your passband does not cover any of the three signal frequency components. This can be seen directly on the graph (the second figure, containing the impulse response, has too fast variations compared with the signal). Or it can be seen from the numbers 0.5760 and 0.3655 in

h = (sin(0.5760*pi*(n-M))-sin(0.3665*pi*(n-M)))./pi./(n-M);

Why did you choose those numbers? The normalized frecuencies of the signal are [2 5 8]/60, that is, 0.0333 0.0833 0.1333. They are all below the passband [.3665 .5760]. As a result, the filter greatly attenuates the three components of the input signal.

Say you want to isolate the center frequency component (f = 5000 Hz, or f/fs = 0.08333 normalized frequency). You need a bandpass filter than lets that frequency through, and rejects the others. So you would use for example a normalized passband [.06 .1]:

h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M);
h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// auto adjustment to avoid the 0/0 sample

A second problem with your code is that it gives two errors because you index h with 0. To solve this, change

n = 0:(N-1);

to

n = 1:N;

and

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j);

to

y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1);

So, the modified code to isolate the central frequency component is:

fs = 60000;
T = 0.006;
t = 0:1/fs:T;
x = 0.3*sin(2*pi*2000*t) + sin(2*pi*5000*t) + 0.4*sin(2*pi*8000*t); 

M = 151;
N = 303;
n = 1:N; %// modified
h = (sin(0.06*pi*(n-M))-sin(0.1*pi*(n-M)))./pi./(n-M); %// modified
 h(n==M) = (h(n==M+1)+h(n==M-1))/2; %// modified


y = zeros(1,length(x)+length(h)-1);
for i = 1:length(x)
for j = 1:length(h)
    y(i+j-1) = y(i+j-1) + x(i)*h(length(h)-j+1); %// modified
end
end

figure(21)
ax1 = subplot(311);
plot(x);
title('Input Signal');
ax2 = subplot(312);
plot(h);
title('FIR');
ax3 = subplot(313);
plot(y);
title('Output Signal');
linkaxes([ax1,ax2,ax3],'x')
ax2.XLim = [0,length(y)];

The result is as follows.

As can be seen, only the central frequency component is present in the output signal.

It's also observed that the envelope of the output signal is not constant. That's because the duration of the input signal is comparable to the filter lenght. That is, you are seeing the transient response of the filter. Note that the envelop rise times is approximately the length of the filter's impulse response h.

It's interesting to note an interesting trade-off here (signal processing is full of trade-offs). To make the transient shorter you could use a wider passband (filter length is inversely proportional to passband). But then the filter would be less selective, that is, it would have less attenuation at the unwanted frequencies. For example, see the result with passband [.04 .12]:

As expected, the transient is now shorter, but the desired frequency is less pure: you can see some "modulation" caused by remains of the other frequencies.