# Thread Subject: Notch filter implementation

 Subject: Notch filter implementation From: vinay Date: 2 Oct, 2010 15:35:04 Message: 1 of 11 Hii im a novice to matlab and im learning it especially the signal processing part I've taken a signal and want to eliminate one frequency from it using notch filter.  I've written some code but the filtered output appears to be same s the input! >> Fs=1000; T=1/Fs; t=(0:Fs-1)*T; x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); subplot(2,1,1); plot(Fs*t(1:50),x(1:50)); wo = 50/(1000/2); bw = wo/35; [b,a] = iirnotch(wo,bw); y=filter(b,a,x); subplot(2,1,2); plot(Fs*t(1:50),y(1:50)); Can anyone help me with this??
 Subject: Notch filter implementation From: Wayne King Date: 2 Oct, 2010 21:37:03 Message: 2 of 11 "vinay " wrote in message ... > Hii im a novice to matlab and im learning it especially the signal processing part > > I've taken a signal and want to eliminate one frequency from it using notch filter. > I've written some code but the filtered output appears to be same s the input! > > >> Fs=1000; > T=1/Fs; > t=(0:Fs-1)*T; > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > subplot(2,1,1); > plot(Fs*t(1:50),x(1:50)); > wo = 50/(1000/2); bw = wo/35; > [b,a] = iirnotch(wo,bw); > y=filter(b,a,x); > subplot(2,1,2); > plot(Fs*t(1:50),y(1:50)); > > Can anyone help me with this?? Hi, I think your notch filter is working pretty well. The problem is that you have a signal with just two sinusoidal components, and the component at 120 Hz is 3 dB larger than the one at 50 Hz. If you look at your filter output in frequency, you will see that the 50 Hz component has been removed nicely. Fs=1000; T=1/Fs; t=(0:Fs-1)*T; x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); subplot(2,1,1); plot(Fs*t(1:50),x(1:50)); wo = 50/(1000/2); bw = wo/60; [b,a] = iirnotch(wo,bw); y=filter(b,a,x); ydft = fft(y); xdft = fft(x); freq = 0:(Fs/length(x)):500; subplot(211) plot(freq,abs(xdft(1:501)).^2) subplot(212) plot(freq,abs(ydft(1:501)).^2); Wayne
 Subject: Notch filter implementation From: vinay Date: 4 Oct, 2010 10:06:08 Message: 3 of 11 "Wayne King" wrote in message ... > "vinay " wrote in message ... > > Hii im a novice to matlab and im learning it especially the signal processing part > > > > I've taken a signal and want to eliminate one frequency from it using notch filter. > > I've written some code but the filtered output appears to be same s the input! > > > > >> Fs=1000; > > T=1/Fs; > > t=(0:Fs-1)*T; > > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > > subplot(2,1,1); > > plot(Fs*t(1:50),x(1:50)); > > wo = 50/(1000/2); bw = wo/35; > > [b,a] = iirnotch(wo,bw); > > y=filter(b,a,x); > > subplot(2,1,2); > > plot(Fs*t(1:50),y(1:50)); > > > > Can anyone help me with this?? > > Hi, I think your notch filter is working pretty well. The problem is that you have a signal with just two sinusoidal components, and the component at 120 Hz is 3 dB larger than the one at 50 Hz. If you look at your filter output in frequency, you will see that the 50 Hz component has been removed nicely. > > Fs=1000; > T=1/Fs; > t=(0:Fs-1)*T; > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > subplot(2,1,1); > plot(Fs*t(1:50),x(1:50)); > wo = 50/(1000/2); bw = wo/60; > [b,a] = iirnotch(wo,bw); > y=filter(b,a,x); > ydft = fft(y); > xdft = fft(x); > freq = 0:(Fs/length(x)):500; > subplot(211) > plot(freq,abs(xdft(1:501)).^2) > subplot(212) > plot(freq,abs(ydft(1:501)).^2); > > Wayne Thank you very much!! But i still don't understand why the two signals (x & y) seem similar when i plot them using plot command..can u help me out
 Subject: Notch filter implementation From: Wayne King Date: 4 Oct, 2010 10:43:22 Message: 4 of 11 "vinay " wrote in message ... > "Wayne King" wrote in message ... > > "vinay " wrote in message ... > > > Hii im a novice to matlab and im learning it especially the signal processing part > > > > > > I've taken a signal and want to eliminate one frequency from it using notch filter. > > > I've written some code but the filtered output appears to be same s the input! > > > > > > >> Fs=1000; > > > T=1/Fs; > > > t=(0:Fs-1)*T; > > > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > > > subplot(2,1,1); > > > plot(Fs*t(1:50),x(1:50)); > > > wo = 50/(1000/2); bw = wo/35; > > > [b,a] = iirnotch(wo,bw); > > > y=filter(b,a,x); > > > subplot(2,1,2); > > > plot(Fs*t(1:50),y(1:50)); > > > > > > Can anyone help me with this?? > > > > Hi, I think your notch filter is working pretty well. The problem is that you have a signal with just two sinusoidal components, and the component at 120 Hz is 3 dB larger than the one at 50 Hz. If you look at your filter output in frequency, you will see that the 50 Hz component has been removed nicely. > > > > Fs=1000; > > T=1/Fs; > > t=(0:Fs-1)*T; > > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > > subplot(2,1,1); > > plot(Fs*t(1:50),x(1:50)); > > wo = 50/(1000/2); bw = wo/60; > > [b,a] = iirnotch(wo,bw); > > y=filter(b,a,x); > > ydft = fft(y); > > xdft = fft(x); > > freq = 0:(Fs/length(x)):500; > > subplot(211) > > plot(freq,abs(xdft(1:501)).^2) > > subplot(212) > > plot(freq,abs(ydft(1:501)).^2); > > > > Wayne > > Thank you very much!! > But i still don't understand why the two signals (x & y) seem similar when i plot them using plot command..can u help me out Hi Vinay, You are just seeing the fact that you have a very simple signal which is dominated by the unattenuated 120 Hz component. That coupled with the fact that you are using a notch filter with a very narrow notch (high Q factor) AND you are using the default -3 dB point as the definition of your bandwidth results in the two signals looking very similar in the time domain. As you can see however, they are different: your notch filter is working. If you reduce the Q factor, i.e. widen the bandwidth, and/or define your bandwidth as a point more than -3 dB down from the peak, you'll start to see more of a difference. Note what happens if you define the bandwidth at -20 dB for the same bandwidth.  Fs=1000;  T=1/Fs;  t=(0:Fs-1)*T;  x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t);  wo = 50/(1000/2); bw = wo/35;  [b,a] = iirnotch(wo,bw,-20);  fvtool(b,a) % view filter magnitude response  y=filter(b,a,x);  plot(x(1:100)); hold on;  plot(y(1:100),'r');   Hope that helps, Wayne
 Subject: Notch filter implementation From: vinay Date: 6 Oct, 2010 04:29:04 Message: 5 of 11 "Wayne King" wrote in message ... > "vinay " wrote in message ... > > "Wayne King" wrote in message ... > > > "vinay " wrote in message ... > > > > Hii im a novice to matlab and im learning it especially the signal processing part > > > > > > > > I've taken a signal and want to eliminate one frequency from it using notch filter. > > > > I've written some code but the filtered output appears to be same s the input! > > > > > > > > >> Fs=1000; > > > > T=1/Fs; > > > > t=(0:Fs-1)*T; > > > > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > > > > subplot(2,1,1); > > > > plot(Fs*t(1:50),x(1:50)); > > > > wo = 50/(1000/2); bw = wo/35; > > > > [b,a] = iirnotch(wo,bw); > > > > y=filter(b,a,x); > > > > subplot(2,1,2); > > > > plot(Fs*t(1:50),y(1:50)); > > > > > > > > Can anyone help me with this?? > > > > > > Hi, I think your notch filter is working pretty well. The problem is that you have a signal with just two sinusoidal components, and the component at 120 Hz is 3 dB larger than the one at 50 Hz. If you look at your filter output in frequency, you will see that the 50 Hz component has been removed nicely. > > > > > > Fs=1000; > > > T=1/Fs; > > > t=(0:Fs-1)*T; > > > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > > > subplot(2,1,1); > > > plot(Fs*t(1:50),x(1:50)); > > > wo = 50/(1000/2); bw = wo/60; > > > [b,a] = iirnotch(wo,bw); > > > y=filter(b,a,x); > > > ydft = fft(y); > > > xdft = fft(x); > > > freq = 0:(Fs/length(x)):500; > > > subplot(211) > > > plot(freq,abs(xdft(1:501)).^2) > > > subplot(212) > > > plot(freq,abs(ydft(1:501)).^2); > > > > > > Wayne > > > > Thank you very much!! > > But i still don't understand why the two signals (x & y) seem similar when i plot them using plot command..can u help me out > > Hi Vinay, You are just seeing the fact that you have a very simple signal which is dominated by the unattenuated 120 Hz component. That coupled with the fact that you are using a notch filter with a very narrow notch (high Q factor) AND you are using the default -3 dB point as the definition of your bandwidth results in the two signals looking very similar in the time domain. As you can see however, they are different: your notch filter is working. > > If you reduce the Q factor, i.e. widen the bandwidth, and/or define your bandwidth as a point more than -3 dB down from the peak, you'll start to see more of a difference. > > Note what happens if you define the bandwidth at -20 dB for the same bandwidth. > > Fs=1000; > T=1/Fs; > t=(0:Fs-1)*T; > x=0.7*sin(2*pi*50*t)+sin(2*pi*120*t); > wo = 50/(1000/2); bw = wo/35; > [b,a] = iirnotch(wo,bw,-20); > fvtool(b,a) % view filter magnitude response > y=filter(b,a,x); > plot(x(1:100)); hold on; > plot(y(1:100),'r'); > > Hope that helps, > Wayne Thank you very much...now i understood where i went wrong Can you advice me where to start to master the signal processing field in matlab?
 Subject: Notch filter implementation From: vinay Date: 19 Jan, 2011 07:25:21 Message: 6 of 11 I've designed a notch filter with different bandwidths and used them on a noise corrupted ecg signal. But the results are not as expected. A signal filtered using a notch filter with bandwidth of 3hz appears to be more refined than the signal obtained using a filter with 1hz bandwidth. Can u please explain..
 Subject: Notch filter implementation From: Wayne King Date: 19 Jan, 2011 11:35:11 Message: 7 of 11 "vinay " wrote in message ... > I've designed a notch filter with different bandwidths and used them on a noise corrupted ecg signal. But the results are not as expected. > > A signal filtered using a notch filter with bandwidth of 3hz appears to be more refined than the signal obtained using a filter with 1hz bandwidth. > > Can u please explain.. I'm afraid that you will need to provide code examples for me to help you. You can simulate a signal. Wayne
 Subject: Notch filter implementation From: vinay Date: 20 Jan, 2011 06:25:20 Message: 8 of 11 Here is the code: I've used MIT-BIH105 record from the physionet database. fs=334; N=length(data);%data is the imported signal i.e rec105 t = ((0:length(data)-1)./fs); t=t'; noise=0.5*sin(2*pi*60*t); input=data+noise; subplot(2,1,1); plot(t,input); title('Noise Corrupted ECG Signal'); xlabel('Time [sec]'); ylabel('Amplitude'); Ys = fft(input)/N; equal_space=linspace(0,.5,N/2); freq = fs*equal_space; Ys = Ys(1:ceil(N)/2); subplot(2,1,2); plot(freq,2*abs(Ys)); title('ECG signal in Frequency Domain'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|'); F_0=60; Delta_F=1; [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b refined1=filter(b,a,input); %filter signal figure; subplot(2,1,1); plot(t,refined); d1=refined1-data; hold;plot(t,data,'r');% red plot is of the pure signal Delta_F=3; [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b refined2=filter(b,a,input); %filter signal d2=refined2-data; subplot(2,1,2); plot(t,refined2); f_iirnotch is a function to find out the coefficients of the notch filter function [b,a] = f_iirnotch (F_0,Delta_F,fs) fs = f_clip (fs,0,fs); F_0 = f_clip (F_0,0,fs/2); Delta_F = f_clip (Delta_F,0,fs/4); r = 1 - (Delta_F*pi/fs);%pole radius theta_0 = 2*pi*F_0/fs; b_0 = abs(1 - 2*r*cos(theta_0) + r^2)/(2*abs(1-cos(theta_0)));%gain factor b = b_0*[1,-2*cos(theta_0),1];%num coefficients a = [1,-2*r*cos(theta_0),r^2];%den coefficients function y = f_clip (x,a,b,k,s) % F_CLIP: Clip a value, or a calling argument, to an interval % % Usage: y = f_clip (x,a,b,k,s) % % Inputs: % x = input to be clipped % a = lower limit of clip range % b = upper limit of clip range % k = an optional integer specifying the number of % the calling argument % s = an optional string specifying the name of the % function being called. % Outputs: % y = x clipped to interval [a,b]. If a <= x <= b, % then y = x. % % See also: F_DEADZONE, F_WAIT, F_PROMPT y = min(x,b); y = max(y,a); if (nargin >= 5) & (y ~= x)     fprintf ('\nCalling argument %d of %s was clipped to [%g,%g].\n',k,s,a,b); end So this is the code..i hope you can figure it out
 Subject: Notch filter implementation From: Wayne King Date: 20 Jan, 2011 13:27:04 Message: 9 of 11 "vinay " wrote in message ... > Here is the code: > I've used MIT-BIH105 record from the physionet database. > > fs=334; > N=length(data);%data is the imported signal i.e rec105 > t = ((0:length(data)-1)./fs); > t=t'; > noise=0.5*sin(2*pi*60*t); > input=data+noise; > subplot(2,1,1); > plot(t,input); > title('Noise Corrupted ECG Signal'); xlabel('Time [sec]'); ylabel('Amplitude'); > Ys = fft(input)/N; > equal_space=linspace(0,.5,N/2); > freq = fs*equal_space; > Ys = Ys(1:ceil(N)/2); > subplot(2,1,2); > plot(freq,2*abs(Ys)); > title('ECG signal in Frequency Domain'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|'); > F_0=60; > Delta_F=1; > [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b > refined1=filter(b,a,input); %filter signal > figure; > subplot(2,1,1); > plot(t,refined); > d1=refined1-data; > hold;plot(t,data,'r');% red plot is of the pure signal > Delta_F=3; > [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b > refined2=filter(b,a,input); %filter signal > d2=refined2-data; > subplot(2,1,2); > plot(t,refined2); > > > f_iirnotch is a function to find out the coefficients of the notch filter > > function [b,a] = f_iirnotch (F_0,Delta_F,fs) > > fs = f_clip (fs,0,fs); > F_0 = f_clip (F_0,0,fs/2); > Delta_F = f_clip (Delta_F,0,fs/4); > r = 1 - (Delta_F*pi/fs);%pole radius > theta_0 = 2*pi*F_0/fs; > b_0 = abs(1 - 2*r*cos(theta_0) + r^2)/(2*abs(1-cos(theta_0)));%gain factor > b = b_0*[1,-2*cos(theta_0),1];%num coefficients > a = [1,-2*r*cos(theta_0),r^2];%den coefficients > > function y = f_clip (x,a,b,k,s) > > % F_CLIP: Clip a value, or a calling argument, to an interval > % > % Usage: y = f_clip (x,a,b,k,s) > % > % Inputs: > % x = input to be clipped > % a = lower limit of clip range > % b = upper limit of clip range > % k = an optional integer specifying the number of > % the calling argument > % s = an optional string specifying the name of the > % function being called. > % Outputs: > % y = x clipped to interval [a,b]. If a <= x <= b, > % then y = x. > % > % See also: F_DEADZONE, F_WAIT, F_PROMPT > > y = min(x,b); > y = max(y,a); > if (nargin >= 5) & (y ~= x) > fprintf ('\nCalling argument %d of %s was clipped to [%g,%g].\n',k,s,a,b); > end > > > So this is the code..i hope you can figure it out Hi Vinay, before I attempt to look at all this code (this was a bit more information than I intended to get :) ), do you have the Filter Design Toolbox? If so, I think you should use iirnotch Wayne
 Subject: Notch filter implementation From: Wayne King Date: 20 Jan, 2011 14:05:08 Message: 10 of 11 "Wayne King" wrote in message ... > "vinay " wrote in message ... > > Here is the code: > > I've used MIT-BIH105 record from the physionet database. > > > > fs=334; > > N=length(data);%data is the imported signal i.e rec105 > > t = ((0:length(data)-1)./fs); > > t=t'; > > noise=0.5*sin(2*pi*60*t); > > input=data+noise; > > subplot(2,1,1); > > plot(t,input); > > title('Noise Corrupted ECG Signal'); xlabel('Time [sec]'); ylabel('Amplitude'); > > Ys = fft(input)/N; > > equal_space=linspace(0,.5,N/2); > > freq = fs*equal_space; > > Ys = Ys(1:ceil(N)/2); > > subplot(2,1,2); > > plot(freq,2*abs(Ys)); > > title('ECG signal in Frequency Domain'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|'); > > F_0=60; > > Delta_F=1; > > [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b > > refined1=filter(b,a,input); %filter signal > > figure; > > subplot(2,1,1); > > plot(t,refined); > > d1=refined1-data; > > hold;plot(t,data,'r');% red plot is of the pure signal > > Delta_F=3; > > [b,a] = f_iirnotch (F_0,Delta_F,fs);%creates coefficient vectors a and b > > refined2=filter(b,a,input); %filter signal > > d2=refined2-data; > > subplot(2,1,2); > > plot(t,refined2); > > > > > > f_iirnotch is a function to find out the coefficients of the notch filter > > > > function [b,a] = f_iirnotch (F_0,Delta_F,fs) > > > > fs = f_clip (fs,0,fs); > > F_0 = f_clip (F_0,0,fs/2); > > Delta_F = f_clip (Delta_F,0,fs/4); > > r = 1 - (Delta_F*pi/fs);%pole radius > > theta_0 = 2*pi*F_0/fs; > > b_0 = abs(1 - 2*r*cos(theta_0) + r^2)/(2*abs(1-cos(theta_0)));%gain factor > > b = b_0*[1,-2*cos(theta_0),1];%num coefficients > > a = [1,-2*r*cos(theta_0),r^2];%den coefficients > > > > function y = f_clip (x,a,b,k,s) > > > > % F_CLIP: Clip a value, or a calling argument, to an interval > > % > > % Usage: y = f_clip (x,a,b,k,s) > > % > > % Inputs: > > % x = input to be clipped > > % a = lower limit of clip range > > % b = upper limit of clip range > > % k = an optional integer specifying the number of > > % the calling argument > > % s = an optional string specifying the name of the > > % function being called. > > % Outputs: > > % y = x clipped to interval [a,b]. If a <= x <= b, > > % then y = x. > > % > > % See also: F_DEADZONE, F_WAIT, F_PROMPT > > > > y = min(x,b); > > y = max(y,a); > > if (nargin >= 5) & (y ~= x) > > fprintf ('\nCalling argument %d of %s was clipped to [%g,%g].\n',k,s,a,b); > > end > > > > > > So this is the code..i hope you can figure it out > > Hi Vinay, before I attempt to look at all this code (this was a bit more information than I intended to get :) ), do you have the Filter Design Toolbox? If so, I think you should use iirnotch > > Wayne Hi Vinay, I see from the earlier posts that you must have iirnotch() and I see that your function is just a reimplementation of the 2nd order IIR notch filter in iirnotch(). The bandwidth of the notch appears to be correct to me. Wayne
 Subject: Notch filter implementation From: vinay Date: 19 Feb, 2011 04:40:04 Message: 11 of 11 > > Hi Vinay, before I attempt to look at all this code (this was a bit more information than I intended to get :) ), do you have the Filter Design Toolbox? If so, I think you should use iirnotch > > > > Wayne > > Hi Vinay, I see from the earlier posts that you must have iirnotch() and I see that your function is just a reimplementation of the 2nd order IIR notch filter in iirnotch(). The bandwidth of the notch appears to be correct to me. > > Wayne Hi Wayne, as u have said i have used the iirnotch() function instead of redefining a new one..Here is the code: fs=334; N=length(data); t = ((0:length(data)-1)./fs); t=t'; input=data+0.5*sin(2*pi*60*t); subplot(2,1,1); plot(t,input); title('Noise Corrupted ECG Signal'); xlabel('Time [sec]'); ylabel('Amplitude'); Ys = fft(input)/N; equal_space=linspace(0,.5,N/2); freq = fs*equal_space; Ys = Ys(1:ceil(N)/2); subplot(2,1,2); plot(freq,2*abs(Ys)); title('ECG signal in Frequency Domain'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|'); w0=(60/(fs/2)); bw=1/(fs/2); [b,a] = iirnotch(w0,bw);%creates coefficient vectors a and b refined1=filter(b,a,input); %filter signal figure subplot(2,1,1) plot(t,refined1) title('Refined ECG signal using filter of bandwidth 1Hz'); subplot(2,1,2) bw=5/(fs/2); [b,a] = iirnotch(w0,bw); refined2=filter(b,a,input); plot(t,refined2) title('Refined ECG signal using filter of bandwidth 5Hz'); My question: Why does the signal filtered using 5hz filter more refined than that using 1hz filter? Thanku Vinay

### Everyone's Tags:

Separated by commas
Ex.: root locus, bode

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.