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');