%f_resp14  Frequency Response, V1.4
%          Displays amplitude ratio 20*log10(AO/AS) in dB and phase
%          in deg versus log10(f), where f is extracted from vS(t)
%          Circuit input (left channel):   vS(t)=AS*cos(2*pi*f*t)
%          Circuit ouptut (right channel): vO(t)=AO*cos(2*pi*f*t+phi)

%          2-27-07, 2-15-07, P. Mathys

fname = input('Enter wav file name: ','s');
[y Fs] = wavread(fname);       %Read recorded wav-file
                               %containing one sweep
T = 0.02;                      %Averaging time constant
vS = y(:,1)';                  %Input vS(t) is left channel
AS = sqrt(2*avgT(vS.^2,Fs,T)); %Amplitude of vS(t)
vScos = vS./AS;                %Normalized cosine from vS(t)
dvScosdt = Fs*[0 diff(vScos)]; %d[vScos]/dt
                               %Derivative of normalized cosine
ffest = sqrt(2*avgT(dvScosdt.^2,Fs,T/2))/(2*pi);
                               %Estimated frequency
vO = y(:,2)';                  %Output vO(t) is right channel
AO = sqrt(2*avgT(vO.^2,Fs,T)); %Amplitude of vO(t)
magHdB = 20*log10(AO./AS);     %Magnitude of H in dB

ord = 1000;                    %Filter order
vSsin = shift90(vScos,Fs,ord); %Shift vScos by -90 deg
a = avgT(2*vO.*vScos,Fs,T);    %a = AO*cos(phi)
b = avgT(2*vO.*vSsin,Fs,T);    %b = -AO*sin(phi)
phi = atan2(-b,a);             %Phase phi of vO(t), modulo 2*pi
phi = unwrap(phi);             %Phase unwrapped (can be left out)

ixdisp = [round(T*Fs):length(ffest)];
                               %Indexes for display
fdisp = [10 10000];            %Frequency range for display
subplot(2,1,1)                 %Selects plot in upper half of
                               %plot window
semilogx(ffest(ixdisp),magHdB(ixdisp),'-b')
                               %Plot |H| in dB versus f
xlim(fdisp)
grid
title('Frequency Response of Second Order LPF (from 20 Hz ... 20 kHz Sweep)')
ylabel('20 log_{10}(|H|) [dB]')

subplot(2,1,2)                 %Selects plot in lower half of
                               %plot window
semilogx(ffest(ixdisp),(180/pi)*phi(ixdisp),'-r')
                               %Plot angle(H) in deg versus f
xlim(fdisp)
grid
ylabel('\angle H [deg]')
xlabel('f [Hz]')

figure(gcf)                    %Brings current figure to foreground

function  y = avgT(x,Fs,T)
%avgT   Compute average value of x(t), using integration window of T seconds
%     >>>>> y = avgT(x,Fs,T) <<<<<
%       where  y     y(t) is integral from t-T to T of x(tau)*dtau
%              x     time function x(t), sampled at rate Fs
%              Fs    sampling rate for x(t) and y(t)
%              T     integration window in seconds

%       2-27-07, 2-13-07, P. Mathys

ixT = round(T*Fs);             %Indexes (or samples) in interval T
xdly = [zeros(1,ixT) x];       %x(t) delayed by T
xdly = xdly(1:length(x));      %Trim to length of x(t)
y = (cumsum(x)-cumsum(xdly))/ixT;
                               %Approx to (1/T)*int_{t-T}^t x(tau)*dtau

function  y = shift90(x,Fs,ord)
%shift90  -90 deg phase shift (Hilbert transform filter)
%       >>>>> y = shift90(x,Fs,ord) <<<<<
%         where
%             y: Filter output y(t), sampling rate Fs
%             x: Filter input x(t), sampling rate Fs
%            Fs: Sampling rate for x(t), y(t)
%           ord: Filter order (must be even integer)

%         2-14-07, P. Mathys

x = reshape(x,1,length(x));    %Make x row vector
ord2 = floor(ord/2);           %Half of filter order
t = [-ord2:ord2]/Fs;           %Time axis for h(t)
h = zeros(size(t));            %Initialize h(t)
ix = find(t~=0);
h(ix) = 1./(pi*t(ix));         %Hilbert transform response
h(ix) = h(ix)-cos(pi*Fs*t(ix))./(pi*t(ix));
                               %H(f) limited to |f|<=Fs/2
y = filter(h,1,[x zeros(1,ord2)])/Fs;  %Filter output y(t)
y = y(ord2+1:end);             %Filter delay compensation