% ### EXaveraging.m ###           01.14.11
% Simple code to demonstrate ideas about spectral vs temporal averaging
% Overview
% 1. Create a sinusiod+noise time waveform
% 2. Re-(under)sample it
% 3. Create time-averaged and spectrally-averaged (i.e., averaging the mag.
% of the Fourier transform)

% stimulus= 0: quantized sinusoid + gaussian noise
% stimulus= 1: % quantized sinusoid (w/ modulating freq.) + gaussian noise

% [NOTE - requires: db.m, rfft.m, irfft.m, and resample.m (if you don't have the signal 
%  processing toolbox, you may need to stick to stimulus= 0)]

clear
% --------------------------------
SR= 44100;         % sample rate {44100}
Npoints= 8192;     % length of fft window (# of points) [should ideally be 2^N] {8192}
f= 3.2;            % 'stimulus' freq. [kHz]
noiseAMP= 1;       % amplitude (re 1) of noise {1}
stimulus= 0;        % different stim. types possible (see above)
repeats= 20;        % # of averages to do
% [slightly more advanced parameters for stimulus= 1]
modDEPTH= 0.0;  % modulation depth in Hz of frequency modulation ** {2}
modFREQ= 5;   % appox. freq. of modulation envelope [Hz] ** {40} % (modulation about center frequency)
% --------------------------------
dt= 1/SR;  % spacing of time steps
nPts= repeats*Npoints;      % total # of points
freq= SR*(0:Npoints/2)./Npoints;    % create a freq. array (for FFT bin labeling)
df = SR/Npoints;    % quantize the 'stim.' freq. (so to have an integral # of cycles)
fQ= ceil(1000*f/df)*df;   % quantized natural freq. [Hz] 
[temp indxF]= max(freq==fQ);    % find freq. array index for sinusoid
t= (0:1/SR:(nPts-1)/SR);    % create an array of time points, Npoints long
% +++
if stimulus==0
    % quantized sinusoid + gaussian noise
    wQ= cos(2*pi*fQ*t) + noiseAMP*randn(size(t,2),1)'; 
elseif stimulus ==1
    % quantized sinusoid (w/ modulating freq.) + gaussian noise
    RND2int=(2*modDEPTH*rand(ceil(modFREQ*max(t)),1))-modDEPTH;
    xx=resample(RND2int,ceil(nPts/ceil(modFREQ*max(t))),1);
    xx=xx(1:size(t,2));
    wQ= cos(t.*(fQ+xx)'*2*pi)+ noiseAMP*randn(size(t,2),1)';    
end
% +++
% plot (entire) time waveform and associated FFT
if 1==1
    figure(1); clf;
    subplot(211); plot(t,wQ);
    hold on; grid on; xlabel('Time [s]'); ylabel('Signal'); title('Entire waveform');
    subplot(212); plot((SR*(0:numel(wQ)/2)./numel(wQ))/1000,db(abs(rfft(wQ))))
    axis([min(freq)/1000 max(freq)/1000 -120 10])
    %xlim([0 t(Npoints)/40])
    grid on; xlabel('Frequency [kHz]'); ylabel('FFT, magnitude [dB]'); title('FFT of entire waveform');
end
% +++
% parse up for time averaging
wAVGtime= zeros(Npoints,1);
for nn=1:repeats
    indx= (nn-1)*Npoints + 1;
    wAVGtime= wAVGtime+ wQ(indx:indx+Npoints-1)';
end
wAVGtime= wAVGtime/repeats;
specT= abs(rfft(wAVGtime)); % mag. spec. for time-averaged waveform
% +++
% parse up for spectral averaging
wAVGspec= zeros(Npoints/2+1,1);
for nn=1:repeats
    clear specTEMP;
    indx= (nn-1)*Npoints + 1;
    specTEMP= abs(rfft(wQ(indx:indx+Npoints-1)'));
    wAVGspec= wAVGspec+ specTEMP;   % Note: this isn't a waveform, but is a complex spectrum!
end
wAVGspec= wAVGspec/repeats;
specS= wAVGspec; % mag. spec. for spectral-averaged waveform
% +++
% plot time-averaged waveform (zoomed in) and FFT
minSPECamp= -70;     % dB min for y-axis of spectra
if 1==1
    figure(2); clf;
    subplot(211); plot(t(1:Npoints),wAVGtime);
    hold on; grid on; axis([0 t(Npoints)/40 -1.5 1.5])
    xlabel('time [s]'); ylabel('signal'); title('Temporally averaged waveform (zoomed-in)')
    subplot(212); plot(freq/1000,db(specT)); hold on; grid on;
    plot(freq(indxF)/1000,db(specT(indxF)),'rx','LineWidth',2)
    axis([min(freq)/1000 max(freq)/1000 minSPECamp 5])
    xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]');
    disp(sprintf('Temporal avg. mag. (at stim. freq.)= %g dB', db(specT(indxF))));
end
% +++
% plot time-averaged waveform and FFT
if 1==1
    figure(3); clf;
    subplot(211); plot(t(1:Npoints),irfft(wAVGspec));   % convert back to time domain
    hold on; grid on; axis([0 t(Npoints)/40 -1.5 1.5])
    xlabel('time [s]'); ylabel('signal'); title('Specrally averaged waveform (zoomed-in)')
    subplot(212); plot(freq/1000,db(specS)); hold on; grid on;
    plot(freq(indxF)/1000,db(specS(indxF)),'rx','LineWidth',2)
    axis([min(freq)/1000 max(freq)/1000 minSPECamp 5])
    xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]');
    disp(sprintf('Spectral avg. mag. (at stim. freq.)= %g dB', db(specS(indxF))));
end