複素フーリエ級数 その３

y(t)=, y(t+2π)=y(t) のフーリエ級数展開。（周期T=2π）
固有周波数f0=1/T=1/2π。

真面目にフーリエ係数を計算すると大変だったのでSymbolic math toolboxにやらせた。
実際、以下のように式変形すれば、フーリエ係数は明らかである。

N=3で完全に元の波を再現できる。

clear; close all; clc;

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

a = 1;
f=@(t) sin(t)+2*cos(3*t)+3*sin(3*t);
y=f(t);
plot(t,y,'DisplayName','original');
hold on;

e = exp(1);
syms x;
x2 = sin(x)+2*cos(3*x)+3*sin(3*x);
c = @(n) 1/T*int(x2*exp((-1j*2*pi*f0*n)*x),-T/2,T/2);
c0 = c(0);

nmaxs=[2,3];
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