Thread Subject:

Notch filter implementation

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??

"vinay " <vinay9023remove.this@gmail.com> wrote in message <i87jf8$8f$1@fred.mathworks.com>...
> 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

"Wayne King" <wmkingty@gmail.com> wrote in message <i888lv$rlm$1@fred.mathworks.com>...
> "vinay " <vinay9023remove.this@gmail.com> wrote in message <i87jf8$8f$1@fred.mathworks.com>...
> > 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

"vinay " <vinay9023remove.this@gmail.com> wrote in message <i8c8ug$m9n$1@fred.mathworks.com>...
> "Wayne King" <wmkingty@gmail.com> wrote in message <i888lv$rlm$1@fred.mathworks.com>...
> > "vinay " <vinay9023remove.this@gmail.com> wrote in message <i87jf8$8f$1@fred.mathworks.com>...
> > > 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

"Wayne King" <wmkingty@gmail.com> wrote in message <i8cb4a$bk3$1@fred.mathworks.com>...
> "vinay " <vinay9023remove.this@gmail.com> wrote in message <i8c8ug$m9n$1@fred.mathworks.com>...
> > "Wayne King" <wmkingty@gmail.com> wrote in message <i888lv$rlm$1@fred.mathworks.com>...
> > > "vinay " <vinay9023remove.this@gmail.com> wrote in message <i87jf8$8f$1@fred.mathworks.com>...
> > > > 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?

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..

"vinay " <vinay9023remove.this@gmail.com> wrote in message <ih63l1$319$1@fred.mathworks.com>...
> 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

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

"vinay " <vinay9023remove.this@gmail.com> wrote in message <ih8kgg$4oj$1@fred.mathworks.com>...
> 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

"Wayne King" <wmkingty@gmail.com> wrote in message <ih9d78$8br$1@fred.mathworks.com>...
> "vinay " <vinay9023remove.this@gmail.com> wrote in message <ih8kgg$4oj$1@fred.mathworks.com>...
> > 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

> > 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