clear all
close all
%
% First, read in a 15 second audio clip.  (Blues Brothers Soul Man)
%
%[Y,FS,NBITS]=wavread('SoulManSnippet.wav');
[Y,FS]=audioread('SoulManSnippet.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)';
%
% axis limits
minYHATS = min(min(20*log10(abs(YHATS))));
maxYHATS = max(max(20*log10(abs(YHATS))));
%
% Plot the specturm.
%
figure(1);
plot(freqs,20*log10(abs(YHATS)));
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
axis([-FS/2, FS/2, minYHATS, maxYHATS]);
title('Spectrum of the original signal')
disp('Figure 1 shows the spectrum of the original signal');
disp(' ')
p=audioplayer(Y,FS);
disp('press enter to continue');
disp(' ')
pause;
disp('Playing original clip');
disp(' ')
play(p)
%disp('press enter to continue');
pause;

%
% Next, filter out every thing above 2Khz.
%
% note abs(freqs) <= 2000 is a boxcar
YHATFILTERED=YHATS(:,1).*(abs(freqs) <= 2000);
YHATFILTERED=YHATS(:,2).*(abs(freqs) <= 2000);
%
% Plot the spectrum of the filtered signal.
%
figure(2);
plot(freqs,20*log10(abs(YHATFILTERED)));
axis([-FS/2, FS/2, minYHATS, maxYHATS]);
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
title('2kHz low pass')
YFILTERED=ifft(ifftshift(YHATFILTERED));
disp('Figure 2 shows the spectrum of the low pass filtered signal');
disp(' ')
p=audioplayer(YFILTERED,FS);
%disp('press enter to continue');
pause;
disp('Playing low pass filtered signal');
disp(' ')
play(p)
%disp('press enter to continue');
pause;

%
% Next, we filter out everything below 1 Khz
%
YHATFILTERED2=YHATS(:,1).*(abs(freqs) >=1000);
YHATFILTERED2=YHATS(:,2).*(abs(freqs) >=1000);
YFILTERED2=ifft(ifftshift(YHATFILTERED2));
%
% Plot the spectrum
%
figure(3);
plot(freqs,20*log10(abs(YHATFILTERED2)));
axis([-FS/2, FS/2, minYHATS, maxYHATS]);
xlabel('frequency Hz');
ylabel('Power (db/Hz)');
title('1 kHz high pass');
disp('figure 3 shows the spectrum of the high pass filtered signal');
disp(' ')
p=audioplayer(YFILTERED2,FS);
%disp('press enter to continue');
pause;
disp('Playing high pass filtered signal');
disp(' ')
play(p)
%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));
p=audioplayer(YFILTERED3,FS);
disp('playing the signal with all 0 phases');
disp(' ')
%disp('press enter to continue');
pause;
play(p)
%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(:,1);
YHATFILTERED4=exp(sqrt(-1)*phase).*YHATS(:,2);
% some distortion got introduced here so remove the complex part
YFILTERED4=real(ifft(ifftshift(YHATFILTERED4)));
p=audioplayer(YFILTERED4,FS);
disp('playing the signal with 15*f phase shift');
disp(' ')
%disp('press enter to continue');
pause;
play(p)
%disp('press enter to continue');
pause;
disp('again')
disp(' ')
disp('signal is shifted by cf radians (cf/2pi cycles) and time length of each cycle is 1/f')
disp('so a 15f phase shift in frequency is a (15f/2pi)*(1/f)=15/2pi = 2.39 second time shift')
disp(' ')
play(p)
%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;
p=audioplayer(YFILTERED5,FS);
%disp('press enter to continue');
disp('playing the signal with 180 degree phase shift');
pause;
play(p)
%disp('press enter to continue');
pause;

