%filter_residual_analysis.m
%FILTER_RESIDUAL_ANALYSIS Compute residual power after filter with various frequencies
%Create a graph of residual vs lowpass filter cutoff frequency.
%D.A. Winter recommends using such a graph to choose filter cutoff freq.
%Tested with wcr0301164.emg(SR=10K) and walkdata20030113rtankleZ.txt(SR=120).
%Uses a Butterworth filter.  W. Rose 20061002, 20080515, 20100419

clear all;

infilename = input('Enter name of file containing the data: ','s');
x=load(infilename);  
N=length(x);
fprintf('Number of points = %i \n',N);
SR = input('Enter sampling rate, in Hertz: ');
p = input('Enter total number of poles (even, <=24) to use for filtering: ');
if p>24 p=24; disp('Setting total poles to 24.'); end;
if mod(p,2)~=0 
    p=p+2-mod(p,2); fprintf('Setting total poles to %d.\n',p);
end
pover2=p/2;                 % single pass filter order
freq=[SR/100:SR/100:SR/2];  % cutoff freqs in Hz
mfreq=freq*2/SR;            % cutoff freqs in "Matlab" frequency units
resid1=zeros(1,length(mfreq));
cutoffcorr=[0.644 0.802 0.863 0.896 0.916 0.929 0.939 0.946 0.952 0.957 0.961 0.964];
%Divide desired cutoff freq by this # to get cutoff freq to use for each 
%stage of the double filter.  Diff #s for diff filter orders from 1 thru 12.
%Values above obtained using C=(sqrt(2)-1)^(1/(2*p)).

i=1;
while mfreq(i)/cutoffcorr(pover2)<=1      %loop through the cutoff freqs
    [b,a]=butter(pover2,mfreq(i)/cutoffcorr(pover2));
    y=filtfilt(b,a,x);
    %for j=1:N              %comment out this calculation of residuals
    %    resid1(i)=resid1(i)+(y(j)-x(j))^2; 
    %end
    %resid1(i)=sqrt(resid1(i)/N);   %this & other way give same results
    resid2(i)=std(y-x,1);
    i=i+1;
end

plot(freq(1:i-1),resid2(1:i-1),'-rx');
titlestr=sprintf('%s   %d pole Butterworth x2',infilename,pover2);
title(titlestr,'Interpreter','none');
xlabel('Low Pass Cut Off Frequency (Hz)');
ylabel('Residual Amplitude'); 
grid on;