% Linear prediction applications
% lp2.m

% Find time-varying filter model for a sound clip



% Load the file
fname='missedit'
%fname='crcmstnc'
%fname='surprise'
%fname='spock2'
%fname='mkmydy'
[x fs bits]=wavread(fname);

% Set the block size
N=300;

% Set the filter order
n=25;

% Find the number of blocks that will fit, and trim input signal
numblocks=floor(length(x)/N);
x=x(1:numblocks*N);

% Reshape the input vector input a matrix
xb=reshape(x,N,numblocks);

% Find the filter coefficients for each block
[a g]=lpc(xb,n);
a = real(a);    % MATLAB v.6 returns complex result 

% Show filter response as a mesh plot
hh=zeros(128,numblocks);
for i=1:numblocks
   [hh(:,i) ff]=freqz(1,a(i,:),128,fs);
end

if 0
mesh(1:numblocks,ff/1000,abs(hh))
axis('ij')
xlabel('block number')
ylabel('frequency (kHz)')
title(['Filter model magnitude response: "' fname '"'])
rotate3d on
disp('press a key...')
pause
end

% Show spectrogram of original image
figure
calspec(x,[],fs)
title('Spectrogram of original sound')
disp('press a key...')
pause


% Show filter response as an image
figure
imagesc(abs(hh)), axis('xy')
title('Spectrogram of filter response')
disp('press a key...')
pause


% Generate an excitation source for the filter
f0=100; ss=blp(40,f0,f0,(length(x)-1)/fs,fs);	% blp (simulates glottal pulses)


% Reshape the input signal
ssb=reshape(ss,N,numblocks);

% Apply the filter
y=zeros(size(ssb));
[y(:,1) z]=filter(g(1),a(1,:),ssb(:,1)); %generate state vector
for k=2:numblocks
   [y(:,k) z]=filter(g(k),a(k,:),ssb(:,k),z);
end

% Reshape the output signal
y=reshape(y,1,N*numblocks);

% Listen to the original vowel sound
soundsc(x,fs)
disp('Original sound...')
pause

% Listen to the synthesized sound
soundsc(y,fs)
disp('Synthesized sound...')

% Show the spectrogram of the synthesized sound
figure
calspec(y,[],fs)
title('Spectrogram of synthesized sound')




