%% Discrete Fourier Transform
%% Fast Fourier Transform (FFT)
% MATLAB computes the discrete Fourier transform (DFT) of a signal using the Fast Fourier Transform (FFT). It can be used to determine the frequency content of a signal.
% Real signals are analog and are described using a continuous function y(t)for t\in (-\infty,\infty). When we measure the signal, we do so for a finite period of time and sample it at discrete times. This gives a digital signal.
% We consider n samples of the signal over a record length of T.

clear
n=256;
T=5;
dt=T/n;
t=(1:n)*dt;
f0=6;
f1=10;
y=sin(2*pi*f0*t)+2*sin(2*pi*f1*t);
plot(t,y)
xlabel('t')
ylabel('y')
title('Signal')
set(gca,'LineWidth',1.2,'FontSize',36)

%
% We seek the frequency content using the discrete Fourier transform. The FFT, or Fast Fourier Transform, was developed to make this computation faster using symmetries within the Fourier transform.

Y=fft(y,n);
n2=n/2;
Power=Y.*conj(Y)/n^2;
f=(0:n2)/(n2*2*dt);
stem(f,2*Power(1:n2+1),'LineWidth',5)
xlabel('Frequency')
ylabel('Power')
title('Periodogram')
set(gca,'LineWidth',1.2,'FontSize',36)

%
%% Spectrogram
% The Fourier transform provides information about the frequency in a signal over the entire time but does not say when a frequency appears. Using a sliding window, one can do this. This is captured in a spectrogram. We consider a chirp signal in which the frequency changes continuously over time.

dt = 0.001;
t = 0:dt:2;
f0 = 50;
f1 = 200;
t1 = 2;
x = chirp(t,f0,t1,f1,'quadratic');

%

sound(10*x,1/dt)

%

plot(t,x)
set(gca,'LineWidth',1.2,'FontSize',36)

%

spectrogram(x,128,120,128,1e3,'yaxis')
colormap jet
set(gca,'LineWidth',1.2,'FontSize',36)
set(gcf,'Position',[2500 100 1500 800])

%
%% Compare to FFT

n=length(t);
xhat=fft(x,n);
PSD = xhat.*conj(xhat)/n;
freq = 1/(dt*n)*(0:n);
L = 1:floor(n/2);

%
% The Fourier transform is the cumulative average of the spectrogram.

figure
plot(freq(L),PSD(L),'LineWidth',2.5)
set(gca,'FontSize',24)
xlabel('Frequency (Hz)')

%

clear
[y,Fs]=audioread('LVBsymph5.mp3');

%
%% Application to Sounds
% Ludwig van Beethoven (1770–1827) Symphony No 5 in C-minor, Op 67, First movement: Allegro con brio from https://archive.org/details/LvB_5_1

T=5;
dt=1/Fs;
t=0:dt:T;
Y=y(Fs:T*Fs);
sound(Y,Fs)

n=length(t);
yhat=fft(Y,n);
PSD = yhat.*conj(yhat)/n;
freq = 1/(dt*n)*(0:n);
L = 1:floor(n/2);
figure
plot(freq(L),PSD(L),'LineWidth',2.5)
set(gca,'FontSize',24)
xlabel('Frequency (Hz)')

%

Nfft=2^9;
spectrogram(Y,Nfft,Nfft-20,Nfft,Fs,'yaxis')
colormap jet
set(gca,'LineWidth',1.2,'FontSize',36)
set(gcf,'Position',[2500 100 1500 800])

%
%% Bird Songs
% Here we look at bird sounds such as those from the sitehttp://earbirding.com/blog/archives/category/spectrograms. The second example is a mocking bird heard in Kure Beach on April 28, 2022.

clear
[y,Fs]=audioread('NorthernCardinal.mp3');
% [y,Fs]=audioread('MockingBird.mp3');
T = length(y)/Fs

%

Tf = 12;
Y=y(1:Tf*Fs);
sound(10*Y,Fs)

%

Nfft=2^8;
spectrogram(Y,Nfft,Nfft-20,Nfft,Fs,'yaxis')
colormap jet
set(gca,'LineWidth',1.2,'FontSize',36)
set(gcf,'Position',[2500 100 1500 800])