Page 144 - MATLAB Recipes for Earth Sciences

138 6 Signal Processing

and

interp1(f,unwrap(phase),1/15) * 15/360
ans =
-5.0000

respectively. Since MATLAB uses causal indexing for filters, the phase

needs to be corrected, similar to the delayed output of the filter. In our
example, we used a filter of the length eleven. We have to correct the

phase by (11-1)/2=5. This suggests a zero phase shift of the filter for both
frequencies.
This also works for recursive filters. Assume a simple sine wave with

period 8 and the previously employed recursive fi lter.
clear
t = (1:100)';
x11 = 2*sin(2*pi*t/8);
b11 = [0.0048 0.0193 0.0289 0.0193 0.0048];
a11 = [1.0000 -2.3695 2.3140 -1.0547 0.1874];

m11 = length(b11);
y11 = filter(b11,a11,x11);

Correct the output for the phase shift introduced by causal indexing and plot
both input and output signals.

y11= y11(1+(m11-1)/2:end-(m11-1)/2,1);
y11(end+1:end+m11-1,1)=zeros(m11-1,1);

plot(t,x11,t,y11)
The magnitude is reduced by

1-max(y11(40:60))/2

ans =
0.6465

which is also supported by the magnitude response
[h,w] = freqz(b11,a11,512);

f= w/(2*pi);
magnitude = abs(h);
plot(f,magnitude)
xlabel('Frequency'), ylabel('Magnitude')