How to account for sample rate in a FIR Digital Filter?

Question

As a learning exercise I am trying to simulate a simple notch filter with Matlab. For test purposes I use Matlab to simulate sinewaves of various frequencies and plot frequency vs. gain so I can see the effect of the “notch”.

I’m using the difference equation $$y(n) = x(n) - 2 \cos\theta\cdot x(n-1) + x(n-2)$$ for this simple notch filter where theta is the notch (say $\pi/3$)

I thought x(n) would be the sample input from the sinewave I generate as test input. I'm just not sure how to account for the sample rate when determining the value of y(n).

From what I expect, and with comparison to using freqz([1 -2*cos(theta) 1],1), I get a good match with a sample rate of 10 and a difference equation as shown below.

$$y(n)= f_s [x(n) -2\cos(\theta/f_s)\cdot x(n-1) + x(n-2)]$$

but, as I change the sample rate there is a huge effect on the gain.

What am I missing? My guess is I just don't understand how to account for a change to the sample rate in the difference equation but I haven't seen a good write up on this anywhere

Can anyone help?

The simple code is shown below. As I change sample rate (from 10 to 100 and 1000) the frequency response appears to be correct but the gain appears to be decreasing by 20dB for each power of 10 I increase the sample frequency. Maybe that’s the clue?

..... CODE .. (You can tell I'm quite new to Matlab also... apologies for style!)

---

fs=10;  %Specify Sampling Frequency

Ts=1/fs; %Sampling period.

fb=100;    % Number of frequencies to try


mb=zeros(fb,2);   %Results   (Max / Freq)

w=0.001;           %Frequency of input wave to be samples  (Don't use DC!)
fb_indx=1;         %Index into the results array


THETA=pi/3;          %For a simple test
while fb_indx < (fb+1)

    Ns=round(2*pi*fs/w);  %Nr of time samples to be plotted.

    t=[0:Ts:Ts*(Ns-1)];  %Make time array that contains Ns elements
                     %t = [0, Ts, 2Ts, 3Ts,..., (Ns-1)Ts]



     mb(fb_indx,1)=w;  %Store current test frequency

     y=zeros(1,Ns);     %Result (y(n))

     x=sin(w*t); %create sampled sinusoids at different frequencies

% Here for the notch

  for n = 3:Ns  % loop for number of samples

     y(n) = fs*(x(n)-(2*cos(THETA/fs)*x(n-1))+x(n-2));  % calculate y

  end

% Obtain max. response
  M=max(y);                %Look for highest peak in result

  mb(fb_indx,2)=20*log10(M);  %and store in 20log10 value

    fb_indx=fb_indx+1;

      w=w+(pi/fb);  %Get next frequency  

  end;

% Plot result


figure % create new figure

subplot(2,2,1)           % first subplot

plot(mb(:,1),mb(:,2));
xlabel('Freq - rad/s');
ylabel('20log10 Magnitude');
axis([0 pi -100 20]);

grid on;

subplot(2,2,2);
plot(x);
title('Input');

subplot(2,2,3); 
plot(y);
title('Output');

Answer 1

The sampling frequency is hidden in your variable $\theta$. The angle $\theta$ determines the frequency of the notch relative to the sampling frequency, i.e.

$$\theta=2\pi\frac{f}{f_s}\tag{1}$$

where $f$ is the notch frequency in Hertz, and $f_s$ is the sampling frequency (in Hertz). So if you want to change the sampling frequency while leaving the notch frequency unchanged, you need to change the angle $\theta$ according to Eq. (1).

The gain of your filter depends on the normalized frequency variable $\theta$. At DC the gain is $2(1-\cos\theta)$, and at Nyquist you get a gain of $2(1+\cos\theta)$. By normalization you can keep one of those two gain factors normalized to $1$, but never both. If you want a gain of $1$ at DC and at Nyquist you need an IIR filter or a higher order FIR filter.