# MATLAB Programming/Advanced Topics/Applications and Examples/Filtering

(Redirected from MATLAB Programming/Filtering)

Filtering is a broad subject. For the MATlab wiki I will focus on how to implement filters. For more on the theory of filtering the reader should reference the Digital Signal Processing wiki book.

## Contents

### The Moving Average Filter

Formula:

${\displaystyle y[n]={\frac {1}{h}}\sum _{p=0}^{h-1}x[n-p]}$

MATLAB implementation(All the code here was intended to be put in an M-file):

 clc;
clear;             % clear all

v=.01
f=100;
fs=5000;
t=0:1/fs:.03
x=sin(2*pi*f*t);                             %original signal

r=sqrt(v)*randn(1,length(t));                %noise
Xw=x+r;                                      %signal plus noise (filter input)
% I have chosen h=3
for n=3:length(Xw),
y(n)=sum(Xw(n-2:n))/3;       %y[n] is the filtered signal
end

plot(y);
hold;

plot(x,'r');                                 %plot the original signal over top the
%filtered signal to see the difference


The moving average filter is simple and effective. One of the things that is a problem is the lag associated with the moving average filter. The more samples used the longer the lag experienced(All filters have lag). How much lag can be tolerated is up to the individual.

### The Kalman Filter

The Kalman filter is a recursive method of combining two estimates to determine the truth. A few parameters that are widely used are the initial conditions or current value and measured data.

Equation:

Example:

 n=100;
sigma=(20/6076);
R=100;
Rm=R+sigma*randn;
Rs(1)=Rm(1);
Cs=sigma^2

for i=2:n
Rm(i)=R+sigma*randn;
alpha=Cs/(Cs+sigma^2);
Rs(i)=Rs(i-1)+alpha*(Rm(i)-Rs(i-1));
Cs=(1-alpha)*Cs;
end


All this code does is take a constant value R and adds noise to it. Then it filters the new signal in an effort to separate the noise from the original signal.

### The discrete Fourier transform

DFT quick reference
What is it? DFT element matlab example and comments
How often do you want to sample? sampling frequency ${\displaystyle f_{s}}$
 fsample=100; sample at 100Hz psample=1/fsample; sample period
For how long do you want to sample? time range ${\displaystyle t=0\ldots t_{\rm {max}}}$
 tmax=10; run from 0 to 10 sec. inclusive
How many samples does that give you? ${\displaystyle N}$
 nsamples=tmax/psample+1; +1 because we include t=0 and t=tmax. nsamples=2^ceil(log2(nsamples)); round up to a power of 2 in length so FFT will work. times=(0:psample:(nsamples-1)*psample)'; Create a column vector of sample times.
How far apart are each of the frequency-domain result points? ${\displaystyle \Delta f={f_{s} \over N}}$
 delf=fsample/nsamples; frequencies=(0:delf:(nsamples-1)*delf)'; Column vector of result frequencies
What signal do you want to sample? input ${\displaystyle x(t)}$
 x=sin(2*pi*10*times)+sin(2*pi*3*times); Make a 10Hz sine wave plus a 3Hz sine wave
What are the results? Fourier transform ${\displaystyle X(f)=\sum _{n=0}^{N-1}x(n)e^{-j2\pi nf/N}}$ fft_x=fft(x, length(x));
What frequencies does the signal have? ${\displaystyle X_{m}(f)=|X(f)|}$ fft_x_mag=abs(fft_x);
What phase relationships? ${\displaystyle X_{\phi }(f)=\arctan \left({\operatorname {Im} (X(f)) \over \operatorname {Re} (X(f))}\right)}$ fft_x_phase=unwrap(angle(fft_x));
How do you view the results?
 plot(frequencies, fft_x_mag); Or, to match the amplitude of the magnitude peak to the amplitude of the sine wave, plot(frequencies, (2/nsamples)*fft_x_mag);
What about the power spectrum? ${\displaystyle X_{p}(f)=|X(f)|^{2}}$
 fft_x_power=fft_x_mag.^2; plot(frequencies, fft_x_power);

### References

Lyons, Richard G. Understanding digital signal processing. Upper Saddle River: Prentice Hall PTR, 2001. ISBN 0-201-63467-8. Chapter 3 discusses the DFT.