function hw7(wavfilename,playoriginal,n)
% EE6820 HW#7
% LPC fit
% and sinwave resynthesis of speech signals
% Usage: hw7('wavefilename', playoriginal (optional,default false), #ofsinusoidstokeep (optional, default 3))
% <xlx@columbia.edu>

if nargin<1
   wavfilename=['d:/courses/audio/wav/mdwh0_sx305.wav'];
   % the "cottage cheese" example
end
if nargin<2
   playoriginal=0; 
end
if nargin<3
   n=3; % # of poles to keep   
end
   
% Data in
[d,rate,nbits]=wavread(wavfilename);
% Downsample to 8 kHz
d = resample(d,8000,rate);
rate = 8000;
if playoriginal == 1
   soundsc(d,rate);
end
% 8th-order LPC modeling
[a,e] = lpcfit(d,8);

% convert to frequencies and magnitudes
freqsyn = zeros(size(a,1),n);
magsyn = zeros(size(a,1),n);

gain = a(:,1);
a1 = [ones(size(a,1),1), a(:,2:size(a,2))];
for i=1:size(a,1)
   poles(i,:)= (roots(a1(i,:)))';
   freqs(i,:) = angle(poles(i,:))*rate/(2*pi);
   mags(i,:) = gain(i)./(1-abs(poles(i,:))); %??
   % find the lowest n poles   
      j=find(freqs(i,:)>0);
      freqtemp=freqs(i,j);
      magtemp =mags(i,j);
      [temp,j]=sort(freqtemp);
   % what if less than n points larger than zero?
   if length(j)>=n
      freqsyn(i,1:n) = freqtemp(j(1:n));
      magsyn(i,1:n) = magtemp(j(1:n));
   else
      freqsyn(i,1:n) = [freqtemp(j),rate*ones(1,n-length(j))/2];
      magsyn(i,1:n) = [magtemp(j),zeros(1,n-length(j))];
   end   
end

sws = synthtrax(freqsyn', magsyn', rate, 128);
figure(1);
subplot(211); specgram(d,128,rate);
subplot(212); specgram(sws,128,rate);
figure(2);
subplot(211); plot((1:length(d))/rate, d);
subplot(212); plot((1:length(sws))/rate,sws);
axis([0,length(sws)/rate,-1,1]);

sws(find(sws>1))=1;
sws(find(sws<-1))=-1;
soundsc(sws,rate);
wavwrite(sws,rate,nbits,'../wav/syntemp.wav');

