%
% First, read in a 20 second audio clip.  (Bruce Cockburn)
%
[Y,FS,NBITS]=wavread('clip.wav');
%
% Next, compute the fft.
%
YHAT=fft(Y);
%
% Shift the fft so it's symmetric around 0 Hz.
%
YHATS=fftshift(YHAT);
fstep=44100/length(Y);
%
% Setup the vector of frequencies.  Note the ' is just to make this a
% column vector like Y, YHAT, etc.
%
freqs=(-22050:fstep:22050-fstep)';
%
% Plot the specturm.
%
figure(1);
plot(freqs,20*log10(abs(YHATS)));
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
print -deps origspect.ps
disp('Figure 1 shows the spectrum of the original signal');
disp('To hear the clip, play clip.wav');
disp('press enter to continue');
pause;
%
% Next, filter out every thing above 2Khz.
%
YHATFILTERED=YHATS.*(abs(freqs) <= 2000);
%
% Plot the spectrum of the filtered signal.
%
figure(2);
plot(freqs,20*log10(abs(YHATFILTERED)));
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
print -deps lowpassspect.ps
YFILTERED=ifft(ifftshift(YHATFILTERED));
wavwrite(YFILTERED,FS,NBITS,'clipb.wav');
disp('Figure 2 shows the spectrum of the low pass filtered signal');
disp('To hear the filtered clip, play clipb.wav');
disp('press enter to continue');
pause;
%
% Next, we filter out everything below 1 Khz
%
YHATFILTERED2=YHATS.*(abs(freqs) >=1000);
YFILTERED2=ifft(ifftshift(YHATFILTERED2));
%
% Plot the spectrum
%
figure(3);
plot(freqs,20*log10(abs(YHATFILTERED2)));
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
print -deps highpassspect.ps
%
% Write out the audio clip.
%
wavwrite(YFILTERED2,FS,NBITS,'clipc.wav');
disp('figure 3 shows the spectrum of the high pass filtered signal');
disp('To hear the high pass filtered clip, play clipc.wav');
disp('press enter to continue');
pause;
%
% Playing with the amplitudes at different frequencies has predictable
% effects, but playing with phase can wreck havoc.  
%
% Here, we set all the phases to 0.  The result isn't anything like the
% original signal!
%
YHATFILTERED3=abs(YHATS);
YFILTERED3=ifft(ifftshift(YHATFILTERED3));
wavwrite(YFILTERED3,FS,NBITS,'clipd.wav');
disp('To hear the version of the signal with all 0 phases, play clipd.wav');
disp('press enter to continue');
pause;
%
% One important way in which phase can be adjusted by a filter is a change
% in phase which is linear with frequency.
%
%
% We will shift each frequency by a phase angle of 15*f.  For example, at
% f=22000 Hz, the phase is shifted by phi=330000 radians, which is 52,521
% cycles, or 2.39 seconds.  Also, at 100 Hz, the phase is shifted by 
% phi=1500 radians, or 238.7 cycles, which is 2.39 seconds.  In fact, the
% signal is shifted in time by 2.39 seconds at every frequency!
%
phase=freqs*15.0;
YHATFILTERED4=exp(sqrt(-1)*phase).*YHATS;
YFILTERED4=ifft(ifftshift(YHATFILTERED4));
wavwrite(YFILTERED4,FS,NBITS,'clipe.wav');
disp('To hear the result of the linear phase adjustment play clipe.wav');
disp('press enter to continue');
pause;
%
% Let's play some more.  What if we shift every phase by pi radians?  This
% is equivalent to multiplying each component of the FFT by exp(i*pi), or
% -1.  Since the FFT and inverse FFT are linear operations,
%       IFFT(-1*FFT(Y))=-Y
%
% Can you detect any difference at all?  
%
YFILTERED5=-Y;
wavwrite(YFILTERED5,FS,NBITS,'clipf.wav');
disp('To hear the 180 degree phase shift, play clipf.wav');