% ### quantizeF3.m ###       04.02.09

% code to show effects/necessity of quantizing freq.


clear; clf;

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

f= 900.0;         % freq. {1000}

SR= 44100;         % sample rate {44100}

Npoints= 32768;     % length of fft window (# of points) [should ideally be 2^N] {8192}

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

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)
df = SR/Npoints;
fQ= ceil(f/df)*df;   % quantized natural freq.

disp(sprintf('specified freq. = %g Hz', f));
disp(sprintf('quantized freq. = %g Hz', fQ));



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

% non-quantized version
%%w= cos(2*pi*f*t);
w= sin(2*pi*f*t);

% quantized version
%%wQ= cos(2*pi*fQ*t);
wQ= sin(2*pi*fQ*t);

% could also apply a window too to the non-quantized version...
wH= hanning(Npoints).*w';



% include Hanning windowed version?
if 1==1
    % plot time waveforms for comparison
    figure(1); clf;
    subplot(211)
    plot(t*1000,w,'o-')
    hold on; grid on;
    plot(t*1000,wQ,'rs-')
    plot(t*1000,wH,'k--d')
    axis([0 5 -1.1 1.1])
    legend('regular vers.','quantized vers.','regular w/ Hanning window')
    xlabel('Time [ms]')
    ylabel('Amplitude')
    title('Note: *regular* v. has arbitrary # of cycles in total FFT window; *quantized* v. has an integer #')
    
    subplot(212)
    plot(t*1000,w,'o-')
    hold on; grid on;
    plot(t*1000,wQ,'rs-')
    plot(t*1000,wH,'k--d')
    axis([t(end-200)*1000 t(end)*1000 -1.1 1.1])
     xlabel('Time [ms]')
    ylabel('Amplitude')
    
    
    % now plot fft of both for comparison
    figure(2); clf;
    plot(freq,db(rfft(w)),'o-','MarkerSize',3)
    hold on;
    plot(freq,db(rfft(wQ)),'rs-','MarkerSize',4)
    %legend('regular version','quantized version')
    xlabel('freq. [Hz]')
    ylabel('magnitude [dB]')
    grid on;
    axis([0 2500 -350 10])
    
    plot(freq,db(rfft(wH)),'k--d','MarkerSize',5)
    legend('regular version','quantized version','regular version w/ Hanning window','Location','NorthEast')
end



% EXCLUDE Hanning windowed version?
if 1==0
    
    % plot time waveforms for comparison
    figure(1); clf;
    plot(t,w,'-')
    hold on; grid on;
    plot(t/1000,wQ,'r--')
    axis([0 5 -1.1 1.1])
    legend('regular version','quantized version')
    xlabel('Time [ms]')
    ylabel('Amplitude')
    
    % now plot fft of both for comparison
    figure(2); clf;
    plot(freq,db(rfft(w)),'o-','MarkerSize',3)
    hold on;
    plot(freq,db(rfft(wQ)),'rs-','MarkerSize',4)
    %legend('regular version','quantized version')
    xlabel('freq. [Hz]')
    ylabel('magnitude [dB]')
    grid on;
    axis([0 2500 -350 10])
    
    legend('regular version','quantized version','Location','East')
    
    
end