% Exploring the convergence of Fourier series
% Using the "Electrolarynx" signal.
%
% Samuli Siltanen Sept 2016

% Parameters for controlling the plots
lwidth = 1;
fsize = 18;
% Choose order of approximation
N = 400;

% Read in the signal
[sgnl0,FS] = audioread('data/glottal_excitation.m4a');
sgnl = sgnl0(140000:145000);

% Take a look
figure(1)
clf
subplot(1,2,1)
plot(sgnl0)
subplot(1,2,2)
plot(sgnl)

% Listen to the signals
%sound(sgnl0,FS)
%sound(sgnl,FS)

% Quadrature points
M = length(sgnl);
h = 2*pi/M;
x = [0:(M-1)]*h;

% Compute the Fourier coefficients
a0      = 1/(2*pi)*h*sum(sgnl);
fseries = a0*ones(size(x));
figure(2)
clf
plot(x,sgnl,'k','linewidth',lwidth)
hold on
plot(x,fseries,'r','linewidth',lwidth)
xlim([0 2*pi])
print -dpng 'images/glottalapprox0.png'
for nnn = 1:N
    coef_c  = 1/pi*h*sum(sgnl.*cos(nnn*x));
    coef_s  = 1/pi*h*sum(sgnl.*sin(nnn*x));
    fseries = fseries +...
        coef_c*cos(nnn*x)+...
        coef_s*sin(nnn*x);
    figure(2)
    clf
    subplot(2,1,1)
    plot(x,sgnl,'k','linewidth',lwidth)
    hold on
    plot(x,fseries,'r','linewidth',lwidth)
    xlim([0 2*pi])
    title(['n=',num2str(nnn)],'fontsize',fsize)
    subplot(2,1,2)
    plot(x,abs(sgnl-fseries),'k','linewidth',lwidth)
    axis([0 2*pi 0 max(abs(sgnl))])
    relerr = norm(abs(sgnl-fseries))/norm(sgnl);
    title(['Relative error=',num2str(round(100*relerr)),'%'],'fontsize',fsize)

    % Write image to file
    fname = ['images/glottalapprox',num2str(nnn),'.png'];
    eval(['print -dpng ',fname])

end
