# Problem with Filtering .... resulting in Aliasing

Started by 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%%%%%%%%%%%%%%%%%%%%%%%%%%

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.
```