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