I do not understand this simple bandstop filter example

Question

I have a periodic signal like this one with period $T=120$:

I would like to apply a bandstop filter to the signal, the stop band should go from 90% of $\frac{1}{T}$ to 110% of $\frac{1}{T}$.

So I run this code in Octave

A=-ones(1,35);
B=ones(1,45);
C=-ones(1,20);
D=ones(1,20);
E=[A B C D];
x=repmat(E,[1 6]);

figure
plot(x,'b');
axis([1 numel(x) -1.5 1.5]);
title('Input signal');
print -dpng input.png


figure
plot(x,'b');
axis([1 numel(x) -1.5 1.5]);


signal_period=numel(E);
signal_frequency=1/signal_period;

sampling_period=1;
sampling_frequency=1/sampling_period;

% "The Nyquist frequency is half the sample rate",
% from http://it.mathworks.com/help/signal/ref/fir1.html#inputarg_Wn
Nyquist_frequency=sampling_frequency/2;

signal_normalized_frequency = signal_frequency/Nyquist_frequency;

b = fir1(40,[0.9*signal_normalized_frequency 1.1*signal_normalized_frequency],'stop');
y = filter(b,1,x);

hold on
plot(y,'r');
title('Input (blue), output (red)');

print -dpng filter_response.png

and I get this filtered signal:

I would expect an attenuation of the signal but I only get a slightly "delayed" signal with some "ringing".

What am I missing?

My ideas:

1. I made a mistake in fir1 with the definition of the lower and upper cutoff frequencies, perhaps I do not understand the normalized frequency concept?
2. Maybe the cutoff frequencies are too low? The frequency response of the filter is the following: and I cannot see the classical "deep" stop band of a notch filter in the magnitude graph.

Answer 1

The problem is that you need the stopband to be very close to 0, which is hard to do; this is why your filter looks like a high-pass. There are two possibilities: reduce the sampling rate, so that 1/T is not so close to 0, or increase the filter order.

I tried increasing the filter order with:

T = 120;
fN = 0.5;
s1 = 0.9/T;
s2 = 1.1/T;
b = fir1(400, [s1, s2]/fN, 'stop');

I don't use Octave any more, but on Matlab I got this filter response with freqz(b):

To try downsampling the signal, you can try commands such as downsample and resample (again, from Matlab, but I seem to recall Octave has equivalents to these). Note that, even in this case, you may need to increase the order to obtain the kind of rejection you require.

Answer 2

fir1 designs your filter using a Hamming window one bin larger than the degree you specify. 41 bins is too small to give you the frequency resolution you want. Try a much larger degree, say 300-600.