複素フーリエ級数その２

y(t)= e^{at}, y(t+2π)=y(t) のフーリエ級数展開。（周期T=2π）
固有周波数f0=1/T。
以下のようにフーリエ級数展開するとき、
y=\displaystyle \sum_{-\infty}^{\infty} c_n e^{2\pi f_0 n t}
係数は以下のように計算される。
c_n=\dfrac{2\cdot (-1)^n \cdot \displaystyle \sinh\left(\dfrac{a}{2}T\right)}{T(a-i(2\pi f_0 n))}

clear; close all; clc;

T=2*pi;
f0=1/T;
dt=0.01;
t=-T/2:dt:T/2;

a = 1;
f=@(t) exp(a*t);
y=f(t);
plot(t,y,'DisplayName','original');
hold on;

e = exp(1);
% c = @(n) sinh(a*pi)/pi*(-1)^n/(a^2+n^2)*(a+1j*n);
c = @(n) 2*(-1)^n*sinh(a/2*T)/T/(a-1j*2*pi*f0*n);
% syms x;
% x2 = exp(a*x);
% c = @(n) 1/T*int(x2*exp((-1j*2*pi*f0*n)*x),-T/2,T/2);
c0 = c(0);

nmaxs=[3,7,10];
for idx=1:length(nmaxs)
    nmax = nmaxs(idx);
    ft=fourier_series_complex(t,f0,c0,c,nmax);
    plot(t,ft,'DisplayName',sprintf('N=%d',nmax));
end
legend();

big;

function ft=fourier_series_complex(t,f0,c0,c,nmax)
ft = 0;
for n=-nmax:nmax
    if n==0
        cn = c0;
    else
        cn = c(n);
    end
    ft = ft + cn*exp(1j*(2*pi*f0*n)*t);
end
end