Page 148 - MATLAB Recipes for Earth Sciences

142 6 Signal Processing

plot(t,xn12,'b-',t,xf12,'r-')
axis([0 200 -4 4])

One might now wish to design a new filter with a more rapid transition from

passband to stopband. Such a fi lter needs a higher order. It needs to have a

larger number of filter weights. We now create a 15-order Butterworth fi lter
as an alternative to the above fi lter.
[b13,a13] = butter(15,0.08/0.5);

[h,w] = freqz(b13,a13,1024);
f = w/(2*pi);

plot(f,abs(h)), grid
xlabel('Frequency')
ylabel('Magnitude')

The frequency response is clearly improved. The entire passband is rela-
tively flat at a value of 1.0, whereas the stopband is approximately zero

everywhere. Next we modify our input signal by introducing a third period
of 5. This signal is then used to illustrate the operation of a Butterworth
bandstop fi lter.
x14 = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
plot(t,x14), axis([0 200 -4 4])
[Pxx,f] = periodogram(x14,[],1024,1);
plot(f,abs(Pxx))

The new Butterworth fi lter is a bandstop fi lter. The stopband of the fi lter is
between the frequencies 0.06 and 0.08. It can therefore be used to suppress
the period of 15 corresponding to a frequency of approximately 0.07.

xn14 = x14 + randn(1,length(t));

[b14,a14] = butter(5,[0.06 0.08]/0.5,'stop');
xf14 = filtfilt(b14,a14,x14);
[Pxx,f] = periodogram(xf14,[],1024,1);

plot(f,abs(Pxx))
plot(t,xn14,'b-',t,xf14,'r-'), axis([0 200 -4 4])

The plots show the effect of this filter. The frequency band between 0.06
and 0.08, and therefore also the frequency of 0.07 was successfully removed
from the signal.