%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Calculate the frequency response function for various filters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  Start by constructing a time series we'd like to evaluate
%    We are designing these so x has unit variance.

ntim = 100;         %  The time series will have 100 points

%  Let's make it AR1:
x = zeros(ntim,1);   %  Initialize x
x(1) = randn(1);     %  Initial condition - random

rho1 = 0.8;           %  One lag autocorrelation
for i = 2:ntim;
  x(i) = rho1*x(i-1)+sqrt(1-rho1^2)*randn(1);
end

%  "Delta" function for the impulse response function:
delt = zeros(1000,1); delt(1) = 1;

%  Filtering:  we're going to try to filter this using several different
%  filters, with similar "cutoff" frequencies.  Note the difference
%  between the cutoff frequency for a butterworth filter (half
%  amplitude), and a running mean (which I'm defining as the first zero
%  crossing).  Let's choose a cutoff of omeg = 0.2 (Nyquist = 0.5)

%  For the rest of this, let's calculate the response functions with 512
%  spectral estimates.
N = 512;


%  Note:  for non-recursive filters - note that in order to make this a
%  centered filter, we define "a" (recursive weights) as zeros, with a 1
%  in the middle.  Type "help freqz" to figure this one out if you're
%  having trouble figuring this out.

%  Choose one of the following filters:

%  Method 1:  Running mean
if 0;
b = ones(1,5)/5;
a = zeros(size(b)); a((length(b)+1)/2) = 1;
filt_title = '5-point running mean';
nonrec = 1;  %  This is a non-recursive filter
end;

%  Method 2:  Succesive 1-2-1 filters
if 0;
b = conv([1 2 1]'/4, [1 2 1]'/4);
a = zeros(size(b)); a((length(b)+1)/2) = 1;
filt_title = '2 Succesive 1-2-1''s';
nonrec = 1;  %  This is a non-recursive filter
end;

%  Method 3:  Lanczos-smoothed filter weights (this is a good
%  non-recursive filter).  We'll use an 11th order filter.  Try playing
%  around with different filter lengths.
if 1;
M = 11;
alph = 0.2;
b = 2*alph*sinc(2*alph*[-M:M]).*sinc([-M:M]/M);
a = zeros(size(b)); a((length(b)+1)/2) = 1;
filt_title = 'Lanczos';
nonrec = 1;  %  This is a non-recursive filter
end;

%  Method 4:  9th order Butterworth filter with cutoff at 0.2
if 0;
[b,a] = butter(9, 2/5);
filt_title = '9th ord. Butterworth (alph = 0.2)';
nonrec = 0;  %  This is NOT a non-recursive filter
%  NOTE!  In this case, the "Power" frequency response function is
%  actually the amplitude response function, because we apply the filter
%  forward, then backwards.
end;


%  Now, get the filtered time series and IRF:
if nonrec
  y1 = filter_nr(b, x);
  irf1 = filter_nr(b, delt);
else
  y1 = filtfilt(b, a, x);
  irf1 = filtfilt(b, a, delt);
end

%  Now, get amplitude response function.
[h, w] = freqz(b, a, N);
f = w/(2*pi);

%  Plot results
figure(1); orient landscape; clf;

%  Plot a subset of the filtered data
ind = 1:100;
subplot(2,1,1);
  plot(ind, x(ind), 'k', ind, y1(ind), 'r');
  grid on
  t1 = title(['Original and Filtered Time Series:  ' filt_title], ...
     'fontsize', 16);  
  
% Plot Impulse Response Function:
ind = 1:20;
subplot(2,2,3);
  plot(ind, irf1(ind), 'k');
  set(gca, 'XLim', [0 ind(end)+1]);
  grid on;
  t1 = title(['IRF:  ' filt_title], 'fontsize', 16);  

%  And, plot the amplitude response function
subplot(2,2,4);
  pp = plot(f, h, 'k', f, conj(h).*h, 'r');
  set(pp(1), 'linewidth', 2);
  grid on;
  axis([0 0.5 -0.4 1.1]);
  grid on;
  t1 = title(['Freq. Resp. Func.:  ' filt_title], 'fontsize', 16);
  ll = legend(pp, 'Amplitude', 'Power');
  set(ll, 'fontsize', 12);