Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Low pass filtering
Date: Sat, 7 Nov 2009 21:25:02 +0000 (UTC)
Organization: The MathWorks Inc
Lines: 74
Message-ID: <hd4oje$coj$1@fred.mathworks.com>
References: <hd4k31$7b7$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1257629102 13075 172.30.248.35 (7 Nov 2009 21:25:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 7 Nov 2009 21:25:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1597503
Xref: news.mathworks.com comp.soft-sys.matlab:583265


"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