MATLAB Central Newsreader - Low pass filtering

Low pass filtering

Dias Papas — Sat, 07 Nov 2009 20:08:02 -0500

I am new in signal processing and I want to filter some time series.
I searched and I found the below code.
It seems to work but I am not sure if it works properly.
And I have two questions.

First: Why the first values of the output 'y' are close to 0? Is it correct? Wouldn't be more logic to be close to the input data 'x'?

And the second one: How can I plot the input data 'X' in the frequency domain with values of Hz, so I will see the frequencies of 'X' and easily decide the cut off frecuency.

If I am in completely wrong thought, please tell me.

The code:

fs = 300;
n = 0:1/fs:0.2;
x = cos(2*pi*30*n) + cos(2*pi*50*n) + sin(2*pi*60*n);
wc = pi*40 / (fs/2);
N = 21;
n = 0:1:N-1;
hd = ideal_lp(wc,N);
w_tet = (boxcar(N))';
h = hd.*w_tet;
y = filter(h,1,x);

fr = 0:0.01:pi;
X = freqz(x,1,fr);
H = freqz(h,1,fr);
Y = freqz(y,1,fr);

figure (1)
subplot(311)
plot(fr,abs(X)); axis([0 3.5 0 40])
subplot(312)
plot(fr,abs(H)); axis([0 3.5 0 2])
subplot(313)
plot(fr,abs(Y)); axis([0 3.5 0 40])

figure (2)
subplot(311)
plot(x)
subplot(312)
plot(h)
subplot(313)
plot(y)

Thanks veeeery much
Dias

Re: Low pass filtering

Wayne King — Sat, 07 Nov 2009 21:25:02 -0500

"Dias Papas" <dias_papas@mathworks.com> wrote in message <hd4k31$7b7$1@fred.mathworks.com>...
> I am new in signal processing and I want to filter some time series.
> I searched and I found the below code.
> It seems to work but I am not sure if it works properly.
> And I have two questions.
>
> First: Why the first values of the output 'y' are close to 0? Is it correct? Wouldn't be more logic to be close to the input data 'x'?
>
> And the second one: How can I plot the input data 'X' in the frequency domain with values of Hz, so I will see the frequencies of 'X' and easily decide the cut off frecuency.
>
> If I am in completely wrong thought, please tell me.
>
> The code:
>
> fs = 300;
> n = 0:1/fs:0.2;
> x = cos(2*pi*30*n) + cos(2*pi*50*n) + sin(2*pi*60*n);
> wc = pi*40 / (fs/2);
> N = 21;
> n = 0:1:N-1;
> hd = ideal_lp(wc,N);
> w_tet = (boxcar(N))';
> h = hd.*w_tet;
> y = filter(h,1,x);
>
> fr = 0:0.01:pi;
> X = freqz(x,1,fr);
> H = freqz(h,1,fr);
> Y = freqz(y,1,fr);
>
> figure (1)
> subplot(311)
> plot(fr,abs(X)); axis([0 3.5 0 40])
> subplot(312)
> plot(fr,abs(H)); axis([0 3.5 0 2])
> subplot(313)
> plot(fr,abs(Y)); axis([0 3.5 0 40])
>
> figure (2)
> subplot(311)
> plot(x)
> subplot(312)
> plot(h)
> subplot(313)
> plot(y)
>
>
> Thanks veeeery much
> Dias

Hi Dias, are you using R2009a or R2009b? If so and you have the Signal Processing Toolbox, how about using fdesign.lowpass? Using your signal, but lengthening it a bit (61 points is a pretty short signal).

n = 0:1/fs:1;
x = cos(2*pi*30*n) + cos(2*pi*50*n) + sin(2*pi*60*n);
H = fdesign.lowpass('Fp,Fst,Ap,Ast',40,45,1,60,300); %passband 40 Hz,stop at 45
% 1 dB passband ripple, 60 dB stopband attenuation, 300 Hz sampling frequency
D = design(H,'equiripple'); % design filter--equiripple FIR design
fvtool(D) %look at filter magnitude response
y = filter(D,x); % Filter input signal

If you look at the help for fdesign.lowpass you will see there are a number of specification strings you can use. A number of them let you specify the filter order. If you want to use the window method for example, you can specify:

H = fdesign.lowpass('N,Fc',50,40,300); %filter order of 50, cutoff of 40 Hz, Fs=300
designmethods(H) % designmethods() tell you what the valid design methods are
% returns window in this case
D = design(H,'window'); %window not necessary here--it's the only method
y = filter(D,x); % filter input signal

In terms of plotting the output iin the frequency domain, can I recommend you read the documentation on fft() and freqz()? If you want to plot the frequency response of your filter, calling freqz() with no output arguments and the proper sampling frequency will produce a plot for you. Alternatively, use two output arguments and plot the absolute value, absolute value squared, or the logarithm of the former as a function of the frequency vector.

To look at the result of filtering your signal, I think you should spend some time reading the documentation for fft(), and or psd(). There a number of worked examples in the Matlab doc.

Hope that helps,
Wayne

Re: Low pass filtering

Dias Papas — Sat, 07 Nov 2009 23:56:02 -0500

"Wayne King" <wmkingty@gmail.com> wrote in message <hd4oje$coj$1@fred.mathworks.com>...
> "Dias Papas" <dias_papas@mathworks.com> wrote in message <hd4k31$7b7$1@fred.mathworks.com>...
> > I am new in signal processing and I want to filter some time series.
> > I searched and I found the below code.
> > It seems to work but I am not sure if it works properly.
> > And I have two questions.
> >
> > First: Why the first values of the output 'y' are close to 0? Is it correct? Wouldn't be more logic to be close to the input data 'x'?
> >
> > And the second one: How can I plot the input data 'X' in the frequency domain with values of Hz, so I will see the frequencies of 'X' and easily decide the cut off frecuency.
> >
> > If I am in completely wrong thought, please tell me.
> >
> > The code:
> >
> > fs = 300;
> > n = 0:1/fs:0.2;
> > x = cos(2*pi*30*n) + cos(2*pi*50*n) + sin(2*pi*60*n);
> > wc = pi*40 / (fs/2);
> > N = 21;
> > n = 0:1:N-1;
> > hd = ideal_lp(wc,N);
> > w_tet = (boxcar(N))';
> > h = hd.*w_tet;
> > y = filter(h,1,x);
> >
> > fr = 0:0.01:pi;
> > X = freqz(x,1,fr);
> > H = freqz(h,1,fr);
> > Y = freqz(y,1,fr);
> >
> > figure (1)
> > subplot(311)
> > plot(fr,abs(X)); axis([0 3.5 0 40])
> > subplot(312)
> > plot(fr,abs(H)); axis([0 3.5 0 2])
> > subplot(313)
> > plot(fr,abs(Y)); axis([0 3.5 0 40])
> >
> > figure (2)
> > subplot(311)
> > plot(x)
> > subplot(312)
> > plot(h)
> > subplot(313)
> > plot(y)
> >
> >
> > Thanks veeeery much
> > Dias
>
> Hi Dias, are you using R2009a or R2009b? If so and you have the Signal Processing Toolbox, how about using fdesign.lowpass? Using your signal, but lengthening it a bit (61 points is a pretty short signal).
>
> n = 0:1/fs:1;
> x = cos(2*pi*30*n) + cos(2*pi*50*n) + sin(2*pi*60*n);
> H = fdesign.lowpass('Fp,Fst,Ap,Ast',40,45,1,60,300); %passband 40 Hz,stop at 45
> % 1 dB passband ripple, 60 dB stopband attenuation, 300 Hz sampling frequency
> D = design(H,'equiripple'); % design filter--equiripple FIR design
> fvtool(D) %look at filter magnitude response
> y = filter(D,x); % Filter input signal
>
> If you look at the help for fdesign.lowpass you will see there are a number of specification strings you can use. A number of them let you specify the filter order. If you want to use the window method for example, you can specify:
>
> H = fdesign.lowpass('N,Fc',50,40,300); %filter order of 50, cutoff of 40 Hz, Fs=300
> designmethods(H) % designmethods() tell you what the valid design methods are
> % returns window in this case
> D = design(H,'window'); %window not necessary here--it's the only method
> y = filter(D,x); % filter input signal
>
> In terms of plotting the output iin the frequency domain, can I recommend you read the documentation on fft() and freqz()? If you want to plot the frequency response of your filter, calling freqz() with no output arguments and the proper sampling frequency will produce a plot for you. Alternatively, use two output arguments and plot the absolute value, absolute value squared, or the logarithm of the former as a function of the frequency vector.
>
> To look at the result of filtering your signal, I think you should spend some time reading the documentation for fft(), and or psd(). There a number of worked examples in the Matlab doc.
>
> Hope that helps,
> Wayne

Wayne really thank you.
Ill study what you ve tell me and if i hope i dont have any further question to bother you!