I acquired the data at 1Khz which is too high for my analysis purposes as i had to see the frequency response down on low frequency range below 25Hz. I have tried using matlab filter topic but when i am using decimate function in matlab (below), i am seeing the aliasing/folding affect. i really want to do just low pass filter the data and downsample to 25hz and take psd/ pwelch to show the trends. i have data of about 44000 points taken at 1khz. Also plz suggest which type of filter is best and how to include it. plz can somebody help me how to perform filtering and down sampling i dont know in details about filters and also a beginner to matlab. plz see .............. the code i wrote is below and help me. Zara %%%%%%%%%code%%%%%%%%%%%%%%%%%%%%%%%%%% FFTData1=load('C:\Documents and Settings\Owner\My Documents\T1.txt') FFTData2=load('C:\Documents and Settings\Owner\My Documents\T2.txt') Fs = 1000; % Sampling frequency T = 1/Fs; % Sample time figure (1) subplot(3,2,1) L1= length(FFTData1) t = (0:L1-1)*T; % Time vector plot(Fs*t(1:L1),FFTData1(1:L1)) xlabel('Time (sec)'); ylabel('Diff'); title('orignal signal T1') grid subplot(3,2,2) L2= length(FFTData2) t = (0:L2-1)*T; % Time vector plot(Fs*t(1:L2),FFTData2(1:L2)); xlabel('Time (sec)'); ylabel('Diff'); title('orignal signal T2'); grid % Detrend/mean is removed from the data ----------------%% X1=FFTData1 X2=FFTData2 m1 = mean(FFTData1) % caculate mean m2 = mean(FFTData2) % caculate mean FFTData1 = DETREND(X1,'m1') FFTData2 = DETREND(X2,'m2') %Detrend/mean Signal--------------------------- Fs = 1000; % Sampling frequency T = 1/Fs; % Sample time subplot(3,2,3) L1= length(FFTData1) t = (0:L1-1)*T; % Time vector plot(Fs*t(1:L1),FFTData1(1:L1),'m') xlabel('Time (sec)'); ylabel('Diff'); title('Detrend/mean removed signal T1') grid subplot(3,2,4) L2= length(FFTData2) t = (0:L2-1)*T; % Time vector plot(Fs*t(1:L2),FFTData2(1:L2),'m') xlabel('Time (sec)'); ylabel('Diff'); title('Detrend/mean is removed T2'); grid %%%%%%%%%%%filtering & sampling%%%%%%%%%%%%%%%%%%%%%%%% d = 20 %with r = Wn = fc/(fs/2) samplingFreq=1000/d; wa1=decimate(FFTData1,15,64,'fir'); wb1=decimate(FFTData2,15,64,'fir'); subplot(3,2,5) plot(wa1,'g') axis([0 length(wa1) -0.1 0.1]) xlabel('Time (sec)'); ylabel('Diff'); title('Filtered signal T1') grid subplot(3,2,6) plot(wb1,'g') axis([0 length(wb1) -0.1 0.1]) xlabel('Time (sec)'); ylabel('Diff'); title('Filtered signal T2'); grid win=512; % Change Window Size lim=win/2; Sampling=samplingFreq; [Pxa,fc] = psd(wa1,win,Sampling,[],50) % Plot the one-sided PSD. fcavg=sum(fc.*Pxa)/sum(Pxa) Pxa = Pxa./Sampling figure(2) subplot(2,1,1) plot(fc,Pxa); xlabel('Frequency (Hz)'); title('T1 by PSD'); ylabel('PSD (-)'); grid [Pxb,fd] = psd(wb1,win,Sampling,[],50) % Plot the one-sided PSD. fdavg=sum(fd.*Pxb)/sum(Pxb) Pxb = Pxb./Sampling subplot(2,1,2) plot(fd,Pxb,'m'); xlabel('Frequency (Hz)'); title('T2 by PSD'); ylabel('PSD (-)'); grid [Pya,fe] = pwelch(wa1,[],[],512,Sampling,'onesided') feavg=sum(fe.*Pya)/sum(Pya) figure(3) subplot(2,1,1) plot(fe,Pya); xlabel('Frequency (Hz)'); title('T1 by pwelch'); ylabel('PSD by pwelch(-)'); grid [Pyb,ff] = pwelch(wb1,[],[],512,Sampling,'onesided') ffavg=sum(ff.*Pyb)/sum(Pyb) subplot(2,1,2) plot(ff,Pyb,'m'); xlabel('Frequency (Hz)'); title('T2 by pwelch'); ylabel('PSD by pwelch (-)'); grid
Problem with Filtering .... resulting in Aliasing
Started by ●November 14, 2008
Reply by ●November 14, 20082008-11-14
I think you are not using the correct value on the decimate function Y = DECIMATE(X,R) resamples the sequence in vector X at 1/R times the original sample rate. The resulting resampled vector Y is R times shorter, i.e., LENGTH(Y) = CEIL(LENGTH(X)/R). By default, DECIMATE filters the data with an 8th order Chebyshev Type I lowpass filter with cutoff frequency .8*(Fs/2)/R, before resampling. Y = DECIMATE(X,R,N) uses an N'th order Chebyshev filter. For N greater than 13, DECIMATE will produce a warning regarding the unreliability of the results. See NOTE below. Y = DECIMATE(X,R,'FIR') uses a 30th order FIR filter generated by FIR1(30,1/R) to filter the data. Y = DECIMATE(X,R,N,'FIR') uses an Nth FIR filter. I think you are trying to decimate by 20. I would recommend decimating by 2 and then by 10. Remember that the R value means the decimation rate. Cheers. On Nov 14, 10:45�am, "zara80" <pinkpeta...@yahoo.co.uk> wrote:> I acquired the data at 1Khz which is too high for my analysis purposes as i > had to see the frequency response down on low frequency range below 25Hz. > > I have tried using matlab filter topic but when i am using decimate > function in matlab (below), i am seeing the aliasing/folding affect. > > i really want to do just low pass filter the data and downsample to 25hz > and take psd/ pwelch to show the trends. i have data of about 44000 points > taken at 1khz. Also plz suggest which type of filter is best and how to > include it. > > plz can somebody help me how to perform filtering and down sampling i dont > know in details about filters and �also a beginner to matlab. > > plz see .............. the code i wrote is below and help me. > > Zara > > %%%%%%%%%code%%%%%%%%%%%%%%%%%%%%%%%%%% > FFTData1=load('C:\Documents and Settings\Owner\My Documents\T1.txt') � > FFTData2=load('C:\Documents and Settings\Owner\My Documents\T2.txt') � > > � � � Fs = 1000; � � � � � � � � � �% Sampling frequency > � � � T = 1/Fs; � � � � � � � � � � % Sample time > > � � � figure (1) > > � � � subplot(3,2,1) > � � � L1= length(FFTData1) > � � � t = (0:L1-1)*T; � � � � � � � �% Time vector > � � � plot(Fs*t(1:L1),FFTData1(1:L1)) > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); > � � � title('orignal signal T1') > � � � grid > > � � � subplot(3,2,2) > � � � L2= length(FFTData2) > � � � t = (0:L2-1)*T; � � � � � � � � � � �% Time vector > � � � plot(Fs*t(1:L2),FFTData2(1:L2)); > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); � � � � > � � � title('orignal signal T2'); > � � � grid � � > > �% Detrend/mean is removed from the data ----------------%% > > � � �X1=FFTData1 > � � �X2=FFTData2 > > � � �m1 = mean(FFTData1) � �% caculate mean � > � � �m2 = mean(FFTData2) � �% caculate mean � > > � � �FFTData1 = DETREND(X1,'m1') � � � � > � � �FFTData2 = DETREND(X2,'m2') � � � � � � � � � > > %Detrend/mean Signal--------------------------- > > � � � Fs = 1000; � � � � � � � � � �% Sampling frequency > � � � T = 1/Fs; � � � � � � � � � � % Sample time > > � � � subplot(3,2,3) > � � � L1= length(FFTData1) > � � � t = (0:L1-1)*T; � � � � � � � �% Time vector > � � � plot(Fs*t(1:L1),FFTData1(1:L1),'m') > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); > � � � title('Detrend/mean removed signal T1') > � � � grid > > � � � subplot(3,2,4) > � � � L2= length(FFTData2) > � � � t = (0:L2-1)*T; � � � � � � � � � � �% Time vector > � � � plot(Fs*t(1:L2),FFTData2(1:L2),'m') > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); � � � � > � � � title('Detrend/mean is removed T2'); > � � � grid � � > > %%%%%%%%%%%filtering & sampling%%%%%%%%%%%%%%%%%%%%%%%% > � � � d = 20 �%with r = Wn = fc/(fs/2) > � � � samplingFreq=1000/d; � > > � � �wa1=decimate(FFTData1,15,64,'fir'); � > � � �wb1=decimate(FFTData2,15,64,'fir'); > > � � � subplot(3,2,5) > � � � plot(wa1,'g') > � � � axis([0 length(wa1) -0.1 0.1]) > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); > � � � title('Filtered signal T1') > � � � grid > > � � � subplot(3,2,6) > � � � plot(wb1,'g') > � � � axis([0 length(wb1) -0.1 0.1]) > � � � xlabel('Time (sec)'); � > � � � ylabel('Diff'); � � � � > � � � title('Filtered signal T2'); > � � � grid � � > > � � � �win=512; � � % Change Window Size � � � � � � � � > � � � �lim=win/2; � � � � � � � � � � � � � � � � � � � � > � � � �Sampling=samplingFreq; > > � � � � [Pxa,fc] = psd(wa1,win,Sampling,[],50) � � % Plot the one-sided > PSD. > � � � � fcavg=sum(fc.*Pxa)/sum(Pxa) > � � � � Pxa = Pxa./Sampling > > � � � � figure(2) � � � �subplot(2,1,1) > � � � � plot(fc,Pxa); > � � � � xlabel('Frequency (Hz)'); � > � � � � title('T1 by PSD'); � > � � � � ylabel('PSD (-)'); � > � � � � grid � � � � > > � � � � [Pxb,fd] = psd(wb1,win,Sampling,[],50) �% Plot the one-sided PSD. > � � � � fdavg=sum(fd.*Pxb)/sum(Pxb) > � � � � Pxb = Pxb./Sampling > > � � � � subplot(2,1,2) > � � � � plot(fd,Pxb,'m'); > � � � � xlabel('Frequency (Hz)'); � > � � � � title('T2 by PSD'); � > � � � � ylabel('PSD (-)'); � > � � � � grid � � � > > � � � � [Pya,fe] = pwelch(wa1,[],[],512,Sampling,'onesided') > � � � �feavg=sum(fe.*Pya)/sum(Pya) > > � � � �figure(3) > > � � � �subplot(2,1,1) > � � � �plot(fe,Pya); > � � � �xlabel('Frequency (Hz)'); � > � � � �title('T1 by pwelch'); � > � � � �ylabel('PSD by pwelch(-)'); � > � � � �grid � � � � > > � � � �[Pyb,ff] = pwelch(wb1,[],[],512,Sampling,'onesided') > � � � �ffavg=sum(ff.*Pyb)/sum(Pyb) > > � � � �subplot(2,1,2) > � � � �plot(ff,Pyb,'m'); > � � � �xlabel('Frequency (Hz)'); � > � � � �title('T2 by pwelch'); � > � � � �ylabel('PSD by pwelch (-)'); � > � � � �grid
Reply by ●November 15, 20082008-11-15