% Evaluates the combined output->input transfer function of a sound card
% linux only
% use a patch cable to connect the output to input

% T.Hilmer UH
%
% 2010.12.02    version 2
%   moved emphasis to analyzing one waveform (swept sine) for various
%   filter orders.
% 2010.12.01    version 1

% Dirac vs Sweep
% Energy is one consideration; dirac over speakers would be impossible.

% Analyzing the combined output->input transfer function of the SBLive, I
% found the results from swept sine to be much more reasonable than dirac.
% That is, the swept sine yielded a very flat input spectrum, empirical
% transfer function, and analytical transfer function for all orders of FIR
% tested (w/ invfreqz).

% TO DO:
%   triangle chirp may create less transients than sawtooth.  Has not been
%   analyzed yet (with other scripts).

% PARAMETRIC MODELING TO DETERMINE HARDWARE FILTER COEFFICIENTS
% IIR: yulewalk
% FIR: lpc
% invfreqz
% has to be called with 'complex'


% LEVELS MUST BE OFF, CAN ALSO HEAR CLICKS.  MUST NEGATE THOSE

version = 3;

addpath ../../fft
clear variables

%------------------------
% Hardware Variables
%------------------------
%'Layla24';
Nch = 2;
dev = alsaSettings('Live',Nch);% ALSA device structure

%------------------------
% Setup Variables
%------------------------
freq = 1E3*[-24 24];% output tone frequency: start, stop
wavestr = 'sawtooth';
wavestr = 'triangle';

% freq = 1E3*[15];% output tone frequency: start, stop
% wavestr = 'tone';
%
% freq = 1E3*[5:5:20];% multitone frequencies
% wavestr = 'multitone';

% !! approximate chirp time.  modified by software to integer # samples
Tc = 3;% approximate chirp time in seconds
Nw = 3;% # number of chirps, defines burst length, suggested 2^N for later FFT

% initialize waveform file:
S = wavefile(dev,wavestr,freq,Tc,Nw);

nonblock = true;% non-blocking execution

period = Tc*(Nw+1);% should be enough time buffer for receive

fname_rx = ['/tmp/lera/' datestr(now,30) '_layla24_rx.bin'];

R = receive(dev,fname_rx,period,nonblock);
% initiate receive before transmit to ensure entire signal is captured
T = transmit(dev,S.fname_tx,period);

Ch = 2; % channel to plot
viewReceive(dev,R.rx_file,60,Ch)
viewReceive(dev,S.fname_tx,60,Ch)

break
%------------------------
% Read Input Data
%------------------------

fid = fopen(R.rx_file,'r');
y = fread(fid,'*int16');
fclose(fid);
y = double(y);
y = complex(y(1:2:end),y(2:2:end)) ./ (2^15-1);


fid = fopen(S.fname_tx,'r');
x = fread(fid,'*int16');
fclose(fid);
x = double(x);
x = complex(x(1:2:end),x(2:2:end)) ./ (2^15-1);

figure
plot(abs(y(1:100000)))% strange behavior with amplitude


%--------------------------------------------------------------------------
% Time-Shift the Input to correspond to Output
%--------------------------------------------------------------------------
[r,lag] = xcorr(x,y,length(S.wav));% no sense lagging more than a period
% [r,lag] = xcov(x,dem,length(wav));% no sense lagging more than a period
% 'coeff','biased',and 'none' vary only by a factor, and decrease w/
% increasing lag (desirable) 'unbiased' is periodic and undesirable.
% For one test, xcorr and xcov game the same lag at peak
[~,j] = max(r);
lag = lag(j);
% phi = angle(r);% slightly more accurate angle later

y = circshift(y,lag);% no phase correction here

if 0 % verify zero time lag
    figure
    plot(abs(x),'k')
    hold on
    plot(abs(y),'b')
end

warning('hack for triangle here')
Ns = length(x)/2;
Fs = dev.Fs;
t = (0:Ns-1)/dev.Fs;% time (s)
% crop excess input (safety buffer)
y = y(1:Ns);
x = x(1:Ns);


f = fftFreq(Ns)*Fs/1000;
X = fft(x);
Y = fft(y);
H = Y./X;

% Y ESSENTIALLY IS THE TRANSFER FUNCTION, SINCE X IS = 1 FOR ALL FREQUENCIES
% PROBABLY DO BETTER JUST AVERAGING Y
% 
% OR PERHAPS AVERAGING H
% 
% spectrogram2(y)
% 
% THIS AVERAGED Y SHOULD WORK WELL WITH INVFREQZ

fr = pi*fftFreq(Ns);

[b,a] = invfreqz(H,fr,'complex',16,16);

h = freqz(b,a,fr);

figure
plot(fftshift(fr)/pi*Fs/1000,fftshift(20*log10(abs(h))),'r')



no = [ 1:50 ];% feedforward orders, i.e. # poles
mo = zeros(size(no));% feedback order, i.e. # zeros
% no zeros is FIR

% preallocate result matrices for analytical filter coefficients
b = zeros(max(no),1);
a = b;
h = nan(Ns,length(no));

for j = 1:length(no)
    disp(j)
    [b(1:no(j)+1,j),a(1:mo(j)+1,j)] = invfreqz(H,fr,'complex',no(j),mo(j));
    h(:,j) = freqz(b(1:no(j)+1,j),a(1:mo(j)+1,j),fr);
end

% started to get unstable (singular matrix) results at filter order = 21 for FIR, started to
% get imaginary components at same order.  Jumps in the coefficient
% magnitudes at 8:9 and 16:17 boundaries.

% I'm starting to think this thing isn't causal.  energy before the
% impulse.... Causality can be recreated by simply lagging the output.
%
% IIR creates phase shift as non-linear function of frequency.  Not
% acceptable for audio.  But careful creation can make linear phase in the
% passband


figure
imagesc(log10(abs(b)))

figure
imagesc(angle(b)),colormap([flipud(jet); jet]),colorbar

break
figure
for j = 1:length(no)
    plot(fftshift(fr)/pi*Fs/1000,fftshift(20*log10(abs(h(:,j)))),'k')
    title(sprintf('Order %.0f',no(j)))
    pause
end


z = filter(b(:,16),1,y);

Xs = Y./h(:,16);

break

figure
plot(fftshift(f),fftshift(20*log10(abs(X))),'.k')
hold on
plot(fftshift(f),fftshift(20*log10(abs(Y))),'.b')
% H(z) is essentially ~Y(z) for the dirac impulse.  H(z) has no visual
% form for swept sine
%     plot(fftshift(f),fftshift(20*log10(abs(H(:,do+1)))),'.r')
plot(fftshift(fr)/pi*Fs/1000,fftshift(20*log10(abs(h))),'r')
legend('X','Y','H','Location','SouthEast')
title([wavestr sprintf(' n=%.0f m=%.0f',n,m)])


% hard to know if these coefficients are near accurate.  Test methods?
% Inverse filter a recording.