Forums

Problem with Filtering .... resulting in Aliasing

Started by zara80 November 14, 2008
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 




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&#2013266080;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 &#2013266080;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') &#2013266080; > FFTData2=load('C:\Documents and Settings\Owner\My Documents\T2.txt') &#2013266080; > > &#2013266080; &#2013266080; &#2013266080; Fs = 1000; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Sampling frequency > &#2013266080; &#2013266080; &#2013266080; T = 1/Fs; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; % Sample time > > &#2013266080; &#2013266080; &#2013266080; figure (1) > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,1) > &#2013266080; &#2013266080; &#2013266080; L1= length(FFTData1) > &#2013266080; &#2013266080; &#2013266080; t = (0:L1-1)*T; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Time vector > &#2013266080; &#2013266080; &#2013266080; plot(Fs*t(1:L1),FFTData1(1:L1)) > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); > &#2013266080; &#2013266080; &#2013266080; title('orignal signal T1') > &#2013266080; &#2013266080; &#2013266080; grid > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,2) > &#2013266080; &#2013266080; &#2013266080; L2= length(FFTData2) > &#2013266080; &#2013266080; &#2013266080; t = (0:L2-1)*T; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Time vector > &#2013266080; &#2013266080; &#2013266080; plot(Fs*t(1:L2),FFTData2(1:L2)); > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080; title('orignal signal T2'); > &#2013266080; &#2013266080; &#2013266080; grid &#2013266080; &#2013266080; > > &#2013266080;% Detrend/mean is removed from the data ----------------%% > > &#2013266080; &#2013266080; &#2013266080;X1=FFTData1 > &#2013266080; &#2013266080; &#2013266080;X2=FFTData2 > > &#2013266080; &#2013266080; &#2013266080;m1 = mean(FFTData1) &#2013266080; &#2013266080;% caculate mean &#2013266080; > &#2013266080; &#2013266080; &#2013266080;m2 = mean(FFTData2) &#2013266080; &#2013266080;% caculate mean &#2013266080; > > &#2013266080; &#2013266080; &#2013266080;FFTData1 = DETREND(X1,'m1') &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080;FFTData2 = DETREND(X2,'m2') &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; > > %Detrend/mean Signal--------------------------- > > &#2013266080; &#2013266080; &#2013266080; Fs = 1000; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Sampling frequency > &#2013266080; &#2013266080; &#2013266080; T = 1/Fs; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; % Sample time > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,3) > &#2013266080; &#2013266080; &#2013266080; L1= length(FFTData1) > &#2013266080; &#2013266080; &#2013266080; t = (0:L1-1)*T; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Time vector > &#2013266080; &#2013266080; &#2013266080; plot(Fs*t(1:L1),FFTData1(1:L1),'m') > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); > &#2013266080; &#2013266080; &#2013266080; title('Detrend/mean removed signal T1') > &#2013266080; &#2013266080; &#2013266080; grid > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,4) > &#2013266080; &#2013266080; &#2013266080; L2= length(FFTData2) > &#2013266080; &#2013266080; &#2013266080; t = (0:L2-1)*T; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080;% Time vector > &#2013266080; &#2013266080; &#2013266080; plot(Fs*t(1:L2),FFTData2(1:L2),'m') > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080; title('Detrend/mean is removed T2'); > &#2013266080; &#2013266080; &#2013266080; grid &#2013266080; &#2013266080; > > %%%%%%%%%%%filtering & sampling%%%%%%%%%%%%%%%%%%%%%%%% > &#2013266080; &#2013266080; &#2013266080; d = 20 &#2013266080;%with r = Wn = fc/(fs/2) > &#2013266080; &#2013266080; &#2013266080; samplingFreq=1000/d; &#2013266080; > > &#2013266080; &#2013266080; &#2013266080;wa1=decimate(FFTData1,15,64,'fir'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080;wb1=decimate(FFTData2,15,64,'fir'); > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,5) > &#2013266080; &#2013266080; &#2013266080; plot(wa1,'g') > &#2013266080; &#2013266080; &#2013266080; axis([0 length(wa1) -0.1 0.1]) > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); > &#2013266080; &#2013266080; &#2013266080; title('Filtered signal T1') > &#2013266080; &#2013266080; &#2013266080; grid > > &#2013266080; &#2013266080; &#2013266080; subplot(3,2,6) > &#2013266080; &#2013266080; &#2013266080; plot(wb1,'g') > &#2013266080; &#2013266080; &#2013266080; axis([0 length(wb1) -0.1 0.1]) > &#2013266080; &#2013266080; &#2013266080; xlabel('Time (sec)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; ylabel('Diff'); &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080; title('Filtered signal T2'); > &#2013266080; &#2013266080; &#2013266080; grid &#2013266080; &#2013266080; > > &#2013266080; &#2013266080; &#2013266080; &#2013266080;win=512; &#2013266080; &#2013266080; % Change Window Size &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;lim=win/2; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;Sampling=samplingFreq; > > &#2013266080; &#2013266080; &#2013266080; &#2013266080; [Pxa,fc] = psd(wa1,win,Sampling,[],50) &#2013266080; &#2013266080; % Plot the one-sided > PSD. > &#2013266080; &#2013266080; &#2013266080; &#2013266080; fcavg=sum(fc.*Pxa)/sum(Pxa) > &#2013266080; &#2013266080; &#2013266080; &#2013266080; Pxa = Pxa./Sampling > > &#2013266080; &#2013266080; &#2013266080; &#2013266080; figure(2) &#2013266080; &#2013266080; &#2013266080; &#2013266080;subplot(2,1,1) > &#2013266080; &#2013266080; &#2013266080; &#2013266080; plot(fc,Pxa); > &#2013266080; &#2013266080; &#2013266080; &#2013266080; xlabel('Frequency (Hz)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; title('T1 by PSD'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; ylabel('PSD (-)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; grid &#2013266080; &#2013266080; &#2013266080; &#2013266080; > > &#2013266080; &#2013266080; &#2013266080; &#2013266080; [Pxb,fd] = psd(wb1,win,Sampling,[],50) &#2013266080;% Plot the one-sided PSD. > &#2013266080; &#2013266080; &#2013266080; &#2013266080; fdavg=sum(fd.*Pxb)/sum(Pxb) > &#2013266080; &#2013266080; &#2013266080; &#2013266080; Pxb = Pxb./Sampling > > &#2013266080; &#2013266080; &#2013266080; &#2013266080; subplot(2,1,2) > &#2013266080; &#2013266080; &#2013266080; &#2013266080; plot(fd,Pxb,'m'); > &#2013266080; &#2013266080; &#2013266080; &#2013266080; xlabel('Frequency (Hz)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; title('T2 by PSD'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; ylabel('PSD (-)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080; grid &#2013266080; &#2013266080; &#2013266080; > > &#2013266080; &#2013266080; &#2013266080; &#2013266080; [Pya,fe] = pwelch(wa1,[],[],512,Sampling,'onesided') > &#2013266080; &#2013266080; &#2013266080; &#2013266080;feavg=sum(fe.*Pya)/sum(Pya) > > &#2013266080; &#2013266080; &#2013266080; &#2013266080;figure(3) > > &#2013266080; &#2013266080; &#2013266080; &#2013266080;subplot(2,1,1) > &#2013266080; &#2013266080; &#2013266080; &#2013266080;plot(fe,Pya); > &#2013266080; &#2013266080; &#2013266080; &#2013266080;xlabel('Frequency (Hz)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;title('T1 by pwelch'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;ylabel('PSD by pwelch(-)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;grid &#2013266080; &#2013266080; &#2013266080; &#2013266080; > > &#2013266080; &#2013266080; &#2013266080; &#2013266080;[Pyb,ff] = pwelch(wb1,[],[],512,Sampling,'onesided') > &#2013266080; &#2013266080; &#2013266080; &#2013266080;ffavg=sum(ff.*Pyb)/sum(Pyb) > > &#2013266080; &#2013266080; &#2013266080; &#2013266080;subplot(2,1,2) > &#2013266080; &#2013266080; &#2013266080; &#2013266080;plot(ff,Pyb,'m'); > &#2013266080; &#2013266080; &#2013266080; &#2013266080;xlabel('Frequency (Hz)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;title('T2 by pwelch'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;ylabel('PSD by pwelch (-)'); &#2013266080; > &#2013266080; &#2013266080; &#2013266080; &#2013266080;grid
Hi 

its me again ........... i tried using the decimate function in steps also
still i am seeing folding effect in signal. 

could any filter expert plz tell exactly how to use decimate correctly or
use any other filter.