% lab 8 Matlab script
% 11.30.2000

pause on;
figurecount=1;

% zero-padding demonstration

  % define time vector
  signallength = 2048;
  n=[0:signallength-1];
  fs=44100;
  Ts=1/fs;

  % define two zero-phase sinusoids closely spaced in frequency
  fhz1=250;
  fhz2=270;
  omega1=fhz1*Ts*2*pi;
  omega2=fhz2*Ts*2*pi;
  x1=0.5.*sin((omega1*n));
  x2=0.5.*sin((omega2*n));

  % create composite signal
  x=x1+x2;
  figure(figurecount);
  subplot(3,1,1);
  plot(x1);
  subplot(3,1,2);
  plot(x2);
  subplot(3,1,3);
  plot(x);
  pause

  % analyze and plot spectrum of composite
  X=fftshift(fft(x)); % negative frequencies shifted left for plotting
  figurecount=figurecount+1;
  figure(figurecount);
  zoomlow=31/64;
  zoomhigh=33/64;
  plot(20.*log10(abs(X(zoomlow*signallength:zoomhigh*signallength)))); 
  title('Spectral magnitude computed from two sinusoids: no zero-padding, zoomed on DC');
  pause

  % create zeropadded version of signal
  zeropadfactor=8;
  xz=[x zeros(1,(signallength*(zeropadfactor-1)))];

  % analyze and plot spectrum of zero-padded version
  XZ=fftshift(fft(xz));
  figurecount=figurecount+1;
  figure(figurecount);
  plot(20.*log10(abs(XZ(zoomlow*signallength*zeropadfactor:zoomhigh*signallength*zeropadfactor)))); 
  title('Spectral magnitude computed from two sinusoids: zero-padded, zoomed on DC');
  pause

  % same signal over a longer period of time
  signallength=2048*8;
  n=[0:signallength-1];
  x1=0.5.*sin((omega1*n));
  x2=0.5.*sin((omega2*n));

  % new composite signal
  x=x1+x2;
  figurecount=figurecount+1;
  figure(figurecount);
  subplot(3,1,1);
  plot(x1);
  subplot(3,1,2);
  plot(x2);
  subplot(3,1,3);
  plot(x);
  pause

  % analyze and plot spectrum of composite
  X=fftshift(fft(x)); % negative frequencies shifted left for plotting
  figurecount=figurecount+1;
  figure(figurecount);
  zoomlow=63/128;
  zoomhigh=65/128;
  plot(20.*log10(abs(X(zoomlow*signallength:zoomhigh*signallength)))); 
  title('Spectral magnitude computed from two sinusoids: longer time basis, zoomed on DC');
  pause

  % create zeropadded version of signal
  zeropadfactor=8;
  xz=[x zeros(1,(signallength*(zeropadfactor-1)))];

  % analyze and plot spectrum of zero-padded version
  XZ=fftshift(fft(xz));
  figurecount=figurecount+1;
  figure(figurecount);
  plot(20.*log10(abs(XZ(zoomlow*signallength*zeropadfactor:zoomhigh*signallength*zeropadfactor)))); 
  title('Spectral magnitude computed from two sinusoids: zero-padded, zoomed on DC');
  pause

% COLA demonstrations

  % triangular window, hop size = 1/2
  figurecount=1;
  w1=[bartlett(64); zeros(32,1)];
  w2=[zeros(32,1); bartlett(64)];
  figure(figurecount);
  subplot(3,1,1);
  stem(w1);
  subplot(3,1,2);
  stem(w2);
  subplot(3,1,3);
  stem(w1+w2);
  pause

  % hamming (raised cosine) window, hop size = 1/2
  figurecount=figurecount+1;
  w1=[hamming(64); zeros(32,1)];
  w2=[zeros(32,1); hamming(64)];
  figure(figurecount);
  subplot(3,1,1);
  stem(w1);
  subplot(3,1,2);
  stem(w2);
  subplot(3,1,3);
  stem(w1+w2);
  pause

% spectrogram demonstrations

  [A,fs,bits]=auread('/usr/ccrma/snd/cburns/small-gong.au');

  % good time resolution
  figurecount=figurecount+1; 
  figure(figurecount);
  subplot(1,1,1);
  specgram(A,512,fs,hamming(512),128);

  % good frequency resolution
  figurecount=figurecount+1;
  figure(figurecount);
  subplot(1,1,1);
  specgram(A,4096,fs,hamming(4096),1024);