
% KAIST CS891 Special Lecture on Music Information Retrieval (MIR)
% Lab1.1: FFT analysis and display
% Kyogu Lee (August 2008)

clear all;
close all;

% Generate a sinusoid
fs = 44100; % sampling rate
amp = 1; % amplitude
freq = 440; % frequency in Hz
T = 1/fs; % sampling period
dur = 2; % duration in seconds
phi = pi/2; % phase in radian

% time-domain signal
x = amp*sin(2*pi*freq*[0:T:dur-T]+phi);

% Fourier transform
nfft = 1024; % fft size
xs = x(1:nfft); % snippet of x
X = fft(xs); % complex numbers
Xdb = 20*log10(abs(X)+0.0000001); % magnitude in dB scale
phs = angle(X); % phase in radians

% display
% time-domain wavform
subplot(311);
ti = 0:T:dur-T; % time index in seconds
plot(ti(1:nfft),xs);
title('x(t)');
xlabel('time (seconds)');
ylabel('amplitude');

% fft magnitude
subplot(312);
fi = 0:fs/nfft:fs-fs/nfft; % freq. resolution = fs/nfft
plot(fi,Xdb,'o-'); xlim([0 fs/2]); % show only up to Nyquist limit (fs/2)
title('Magnitude response');
xlabel('frequency (Hz)');
ylabel('magnitude (dB)');

% phase angle
subplot(313);
plot(fi,phs), xlim([0 fs/2]);
title('Phase response');
xlabel('frequency (Hz)');
ylabel('phase (rad)');