clc;
clear all;
close all;
fclose all;


ref=wavread('swisha.wav'); % this is the original full-bandwidth probe tone.
% Do not sample rate convert, or otherwise play with it.

signal=wavread('foo.wav');  %this reads a stereo file, channel 1 i.e. (:,1)
% is the captured signal, channel 2 is the loopback signal 

len=2*length(ref) % set the analysis length. This is good for anything that has a T60 under 5 or 6 seconds.

slen=length(signal) % how long is the signal length. It should be about "len" long.

reft=fft(ref,len);  % transform of probe tone, or reference.

lb=signal(1:slen,2); % loopback signal
imp=signal(1:slen,1); % acoustic test signal from mic.

lbt=fft(lb,len);  %transform loopback at length "len"  n.b. not sure what happens if slen is too long!

impt=fft(imp,len); % transform of actual accosutic mic signal

lbtd=lbt ./ reft; %deconvolve spectrum of loopback

imptd=impt ./ reft; % and deconvolve spectrum of mic signal

close all;  % clear up anything hiding about the graphics system

lbtddb=abs(lbtd(1:len/2+1)); %loop back amplitude spectrum
t=max(lbtddb) % normalize to zero
lbtddb=lbtddb/t; % finish normalization
lbtddb=max(lbtddb, .001); % no fair taking log zero
lbtddb=20*log10(lbtddb); % convert to dB
f=(0:len/2)*22050/(len/2); % calculate frequency scale
f=f'; %make octave happy with it.
subplot(3,1,1) % first plot setup
plot(f,lbtddb); % plot frequency response of loopback signal.
%That ought to be very, very good

subplot(3,1,2);  % set up second plot
imptdb=abs(imptd(1:len/2+1)); % mic capture amplitude spectrum
t=max(imptdb); % get max
imptdb =imptdb /t;%finish normalization
imptdb=max(imptdb, .001);% no log zer0
imptdb=20*log10(imptdb); % convert to dB
plot(f,imptdb); % and plot frequency response of actual system.

%you could use a log scale for frequency if you limit it to 20-22050Hz.
%you could also smooth the frequency response, either via fixed or part-octave

lbat(1:len)=0;  %set up analytic envelope for loopback
lbat(1:len/2+1)=lbtd(1:len/2+1); %load positive frequency part
lba=abs(ifft(lbat)); %abs of ifft gives analytic envelope
t=max(lba); %get max tonormalize
lba=lba/t; %finish norm
lba=max(lba,.01);%limit dynamic range
lba=log10(lba); %take log to see details

subplot(3,1,3)  %set up third plot


sat(1:len)=0; %calculate analytic envelope, log, etc, for mic signal
sat(1:len/2+1)=imptd(1:len/2+1);
sa=abs(ifft(sat));
t=max(sa);
sa=sa/t;
sa=max(sa,.01);
sa=log10(sa);


[junk1,inds]=max(sa);  % inds finds max value for signal
[junk,indx]=max(lba); % and indx max value location for loopback
xx=-100:1000; %set up horizontal scale
xx=xx/44100*1153; % convert feet
plot(xx, sa(indx-100:indx+1000),'g'); %plot  mic signal envelope
hold on; % don't erase the plot
plot(xx, lba(indx-100:indx+1000),'r'); %plot loopback signal envelope

sdif=inds-indx
dist=sdif/44100*1153 % time delay at 72F in feet.