% ### SOAEplayback.m ###    C. Bergevin (7.14.09)

% simple code to play back audio of (filtered) SOAE time waveforms

clear

% -----------------------------------------------
% specifiable parameters
file='./RRrearSOAEwf1B.txt';  % SOAE time waveform (assumes SR=22050 Hz)

%nPts=8192; % # of points in FFT buffer (default=8192)
SR=22050;  % sample rate (needed to determine freq. bins)

% Filter Type
filt= 2;   % high-pass filter (1) or band-pass filter

% High-Pass Filter
cutoffF= 1700;  % cutoff freq. 

% Band-Pass Filter (Recursive Exponential)
CF= 2206;   % filter center freq.
BW= 5;     % ~bandwidth of filter

gain= 40;  % dB boost to signal playback 

% -----------------------------------------------

wf=load(file);

M= size(wf,1);  % # of points in time waveform

% create a time array (if needed)
time=[0:1/SR:M/SR];
time=time(1:end-1);
disp(sprintf('file played back: %s', file));
disp(sprintf('duration of time waveform = %g s', M/SR));

% determine freq. bins for spectrum
freq= [0:SR/M:SR/2];


Sfor=rfft(wf);     % convert to fourier domain (keep complex!)


% apply filter
if filt==1
    temp= find(freq>=cutoffF);
    N1= temp(1);
    % tamp down values outside (hard) window
    Sfor(1:N1)=10^-7;
    disp(sprintf('high pass filter: fc=%g Hz', cutoffF));
elseif filt==2
    % determine filter from specified parameters
    F= recexp2(freq,CF,BW,10);
    Sfor= F.*Sfor';   % apply filter in spectral domain
    disp(sprintf('band-pass filter: CF=%g Hz, BW=%g Hz', CF,BW));
end

% visualize to ensure things look okay
figure(2)
clf
plot(freq,db(rfft(wf)))
hold on;
plot(freq,db(Sfor),'r')
axis([0 7000 -150 -50])
xlabel('freq. [Hz]')
ylabel('mag. [dB arb.]');

% now convert back to a time signal
HPfiltSPEC= irfft(Sfor);



disp(sprintf('starting signal playback'));
% now play back to signal
sound(HPfiltSPEC*(10^(gain/20)),SR,24)
disp(sprintf('playback done'));
