% ### specREP2.m ###       09.17.13


% code to just fiddle with basics of the FFT and what the spectrum of
% various signals look like [updated for BPHS 4090 F2013]

% *****************************************
% Some questions to think about....

% QUESTION: What is fft doing that rfft is not? What sort of symmetry is
% going on? Is there some information that is 'extraneous'? What does it
% mean to have 'negative frequencies'?

% QUESTION: What is the 'dB' scale for the vertical axis of the magnitude?
% Why use a log scale?

% QUESTION: For a click, how does the phase of the unwrapped response
% relate to the time waveform? How does the spectrum change with changing
% the 'duration' of the click (i.e., CLKoff-CLKon)? Forexample, why do
% notches start to manifest? ANd what happens when there a 'a lot' of
% notches?

% QUESTION: What changes can you make such that only the magnitude is
% affected? Only the phase?

% QUESTION: What effects does changing the SR have? # of point in the FFT
% window?

% QUESTION: Why are there differences between rfft and fft output? What is
% the correct way to handle things?

% QUESTION: What does it mean to 'quantize' the frequencies?

% QUESTION: What is different between the different types of noise?


% ------
% Stimulus Type Legend
% stimT= 0 - non-quantized sinusoid
% stimT= 1 - quantized sinusoid
% stimT= 2 - one quantized sinusoid, one un-quantized sinusoid
% stimT= 3 - two quantized sinusoids
% stimT= 4 - click I.e., an impulse)
% stimT= 5 - noise (flat mag.)
% stimT= 6 - chirp (flat mag.)
% stimT= 7 - noise (Gaussian; flat spectrum, random phase)
% stimT= 8 - exponentially decaying sinusoid

clear; clf;

% --------------------------------
% User Inputs

SR= 44100;         % sample rate [Hz]

Npoints= 8192;     % length of fft window (# of points) [should ideally be 2^N]
% [time window will be the same length]

% Stimulus Type (see legend above)
stimT= 0;

f= 1500.0;         % Frequency (for waveforms w/ tones) [Hz]

% Note: Other stimulus parameters (e.g., click properties, two tones) have
% parameters that can be changed below

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

dt= 1/SR;  % spacing of time steps

% create a freq. array (for FFT bin labeling)
freq= [0:Npoints/2];
freq= SR*freq./Npoints;

% quantize the freq. (so to have an integral # of cycles)
% [see quantizeF.m for further aspects on this]
df = SR/Npoints;
fQ= ceil(f/df)*df;   % quantized natural freq.

% create an array of time points, Npoints long
%t= linspace(0,Npoints/SR,Npoints);
t=[0:1/SR:(Npoints-1)/SR];

% compute stimulus
if stimT==0
    % non-quantized sinusoid
    signal= cos(2*pi*f*t);
    disp(sprintf(' \n *Stimulus* - (non-quantized) sinusoid, f = %g Hz \n', f));
    disp(sprintf('specified freq. = %g Hz', f));
    disp(sprintf('quantized freq. = %g Hz', fQ));
    
elseif stimT==1
    % quantized sinusoid
    signal= cos(2*pi*fQ*t);
    disp(sprintf(' \n *Stimulus* - quantized sinusoid, f = %g Hz \n', fQ));
    disp(sprintf('specified freq. = %g Hz', f));
    disp(sprintf('quantized freq. = %g Hz', fQ));
    
elseif stimT==2
    % one quantized sinusoid, one un-quantized sinusoid
    ratio= 1.22;    % specify f2/f2 ratio
    signal= cos(2*pi*fQ*t) + cos(2*pi*ratio*fQ*t);
    disp(sprintf(' \n *Stimulus* - two sinusoids (one quantized, one not) \n'));
elseif stimT==3
    % two quantized sinusoids
    ratio= 1.22;    % specify f2/f2 ratio
    fQ2= ceil(ratio*f/df)*df;
    signal= cos(2*pi*fQ*t) + cos(2*pi*fQ2*t);
    disp(sprintf(' \n *Stimulus* -  two sinusoids (both quantized) \n'));
elseif stimT==4
    % click
    CLKon= 1000;     % index at which click turns 'on' (starts at 1)
    CLKoff= 1001;   % index at which click turns 'off'
    clktemp1= zeros(1,Npoints);
    clktemp2= ones(1,CLKoff-CLKon);
    signal= [clktemp1(1:CLKon-1) clktemp2 clktemp1(CLKoff:end)];
    disp(sprintf(' \n *Stimulus* - Click \n'));
elseif stimT==5
    % noise (flat)
    signal= rand(1,Npoints);
    disp(sprintf(' \n *Stimulus* - Noise1 \n'));
elseif stimT==6
    % chirp (flat)
    f1S= 2000.0;     % if a chirp (stimT=2) starting freq. [Hz] [freq. swept linearly w/ time]
    f1E= 4000.0;   % ending freq.  (energy usually extends twice this far out)
    f1SQ= ceil(f1S/df)*df;      %quantize the start/end freqs. (necessary?)
    f1EQ= ceil(f1E/df)*df;
    % LINEAR sweep rate
    fSWP= f1SQ + (f1EQ-f1SQ)*(SR/Npoints)*t;
    % quantize the freqs. [NOTE: this seems to have a big effect!!]
    %         for mm=1:nPts
    %            fSWP(mm)= ceil(fSWP(mm)/df)*df;
    %         end
    %%%signal = ChirpVolt*sin(2*pi*fSWP.*timeW)';
    signal = sin(2*pi*fSWP.*t)';
    disp(sprintf(' \n *Stimulus* - Chirp \n'));
elseif stimT==7
    % noise (Gaussian)
    Asize=Npoints/2 +1;
    % create array of complex numbers w/ random phase and unit magnitude
    for n=1:Asize
        theta= rand*2*pi;
        N2(n)= exp(i*theta);
    end
    N2=N2';
    % now take the inverse FFT of that using Chris' irfft.m code
    tNoise=irfft(N2);
    % scale it down so #s are between -1 and 1 (i.e. normalize)
    if (abs(min(tNoise)) > max(tNoise))
        tNoise= tNoise/abs(min(tNoise));
    else
        tNoise= tNoise/max(tNoise);
    end
    signal= tNoise;
    disp(sprintf('\n *Stimulus* - Noise2 \n'));
elseif stimT==8
    % exponentially decaying cos
    alpha= 500;
    %signal= exp(-alpha*t).*cos(2*pi*fQ*t);
    %signal= exp(-alpha*t).*cos(2*pi*f*t);
    signal= exp(-alpha*t).*sin(2*pi*fQ*t);
end



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

% *******
% plot time waveform of signal
figure(1); clf
plot(t*1000,signal,'k.-','MarkerSize',5)
grid on; hold on;
xlabel('Time [ms]')
ylabel('Signal')
title('Time Waveform')



% *******
% now plot rfft of the signal

sigSPEC= rfft(signal);

% MAGNITUDE
figure(2); clf;
subplot(211)
plot(freq/1000,db(sigSPEC),'ko-','MarkerSize',3)
hold on; grid on;
ylabel('Magnitude [dB]')
title('Spectrum')
% PHASE
subplot(212)
% radians (wrapped)
%plot(freq,angle(rfft(signal)),'o-','MarkerSize',3)
% cycles (unwrapped)
plot(freq/1000,cycs(sigSPEC),'ko-','MarkerSize',3)
xlabel('Frequency [kHz]')
ylabel('Phase [cycles]')
grid on;
%axis([0 2500 -350 10])




% *******
% now plot fft instead (i.e., include negative frequencies); note that the 'frequency' needs to be
% different for this case (i.e., 'negative' frequencies); phase is not
% unwrapped
if 1==1
    % kludge; correct way to specify?
    %freq2= [-flipud(freq(3:end)')' freq];
    
    NFFT = 2^nextpow2(Npoints);
    %freq2 = SR/2*linspace(0,1,NFFT/2+1);
    freq2 = SR/2*linspace(-1,1,NFFT);
    
    sigSPEC2= fft(signal);
    
    % MAGNITUDE
    figure(3); clf;
    subplot(211)
    plot(freq2,db(abs(sigSPEC2)),'ko-','MarkerSize',3)
    hold on; grid on;
    ylabel('magnitude [dB]')
    % PHASE
    subplot(212)
    % radians (wrapped)
    plot(freq2,angle(sigSPEC2),'ko-','MarkerSize',3)
    % cycles (unwrapped)
    %plot(freq,cycs(rfft(signal)),'o-','MarkerSize',3)
    xlabel('Frequency?? [Hz]')
    ylabel('phase')
    grid on;
    %axis([0 2500 -350 10])
end

% -------
% play the stimuli
sound(signal,SR)


% -------
% compute inverse Fourier transform and plot
if 1==1
    figure(1);
    signalINV= irfft(sigSPEC);
    plot(t*1000,signalINV,'rx','MarkerSize',4)
end