%
%  [amplitude,phase]=filtfiltfreqz(B,A,f)
%
%  Computes and plots the frequency response of a filtfilt filter.  
%
%  Inputs:
%     A,B       filter coefficients.
%     f         Values of the normalized frequency (f/fs) at which
%               to evaluate the response.
%
%  Outputs:
%     amplitude Amplitude of the response at different frequencies.
%
function amplitude=filtfiltfreqz(B,A,f)
%
% First, compute the equivalent s and z values for f.
%
s=2*pi*i*f;
zinv=exp(-s);
%
% We'll need to flip the order of the coefficients in A and B to
% use polyval, since it has the opposite convention from filtfilt.
%
Aflipped=fliplr(A);
Bflipped=fliplr(B);
%
% Now, for each frequency, evaluate Phi(z).
%
Phi=polyval(Bflipped,zinv)./polyval(Aflipped,zinv);
%
% Plot the amplitude of the filtfilt filter.
%
plot(f,20*log10(abs(Phi).^2),'k-');
hold on
xlabel('f/fs');
ylabel('abs(Phi(f))^2 (dB)');
%
% Now, return abs(Phi)^2.
%
amplitude=abs(Phi).^2;

