clear;
fs=40000; %sampling rate
tt=5; %signal duration in seconds
f1=20; %starting frequency
f2=20000; %ending frequency
m=(f2-f1)/tt; %rate of sweep
psi=0; %starting phase
t=0:1/fs:tt; n=0:fs*tt;
%x=zeros(1,tt*fs+1); %initialize signal
x= cos(pi*m*t.*t + 2*pi*f1*t + psi); %synthesize chirp
y=cos(2*pi*f1*t);
subplot(111), 
%stem(n,x,'y'); axis tight;
%hold on; 
subplot(211), plot(n,x); 
xlabel('discrete-time, n'); ylabel('x(n)'); 
title(['Chirp signal, m= ',num2str(m),' kHz, f1= ',num2str(f1),' kHz']);
z=x+y;
wavwrite(0.9*z/(max(abs(z))),fs, 'chirp');
soundsc(z,fs)
Nw=round(length(x)/6);
xw=x(round(length(x)/2)-round(Nw/2): round(length(x)/2)+round(Nw/2));
N=1024;
fx=fft(xw,N); fx=fx(1:N/2); fx=abs(fx);
f=linspace(0,fs/2,N/2);
subplot(212), plot(f,fx); 
figure(gcf)


% Equal Loudness in Time Domain

curve = iso226(60);  %29 points of fletcher munson curve
c = interp(curve, 6896);    %expand 29 point to 200000 points
c = [c zeros(1,17)];

ELchirp = c.*z;    %multiply chirp by equal loudness curve
soundsc(ELchirp, fs);

figure 
subplot(211), plot(curve);
axis tight;
xlabel('time (samples)');  ylabel('Loudness (phons)');
title('Equal Loudness Curve');

subplot(212), plot(ELchirp);
axis tight;
xlabel('time (samples)');  ylabel('Amplitude');
title('Equal Loudness Chirp');

wavwrite(0.9*ELchirp/(max(abs(ELchirp))),fs, 'eqloudchirp');

%  Equal Loudness in Frequency Domain

curve = iso226(60);  %29 points of fletcher munson curve

[vega, fs] = wavread('vega.wav');
soundsc(vega, fs);
L = length(vega);
nfft = 2^nextpow2(L);
f = fs/2*linspace(0,1,nfft/2);
vega = vega';

cc = interp(curve, floor(nfft/29));       %expand 29 points to 512 points 
cc = [cc zeros(1, nfft-length(cc))];

vega=[vega zeros(1,length(vega))];  %append zeros at beginning and end of frame

v = fft(vega, nfft);     % fft of vega.wav
v = v(1:length(v)/2);

ELclip = cc.*v;          %filter fft w/ curve
ELvega = ifft(ELclip, L);

soundsc(ELvega,fs);

figure 
subplot(311), plot(f,2*abs(v(1:nfft/2+1)));
axis tight;
xlabel('frequncy (bins)');  ylabel('Amplitude');
title('Vega Fourier Transform');


subplot(312), plot(vega);
axis tight;
xlabel('time (samples)');  ylabel('Amplitude');
title('Vega Sound Clip');

index = 1:length(ELvega);
subplot(313), plot(index, ELvega);
axis tight;
xlabel('time (samples)');  ylabel('Amplitude');
title('Equal Loudness Sound Clip');

wavwrite(0.9*ELvega/(max(abs(ELvega))),fs, 'eqloudvega');







