* ********************************************************************
* Symbolic Fourier Transform in Matlab (Octave)
*
* f(t) = (1/(2*PI)) * integral(F(w) e^(jwt) dw)
*
* F(w) =              integral(f(t) e^(-jwt) dt)
* ********************************************************************

pkg load symbolic

syms t

PI = sym(pi)
ONE = sym(1)
ENM = sym(e)
JAY = sym(j)

g(t) = sin(4 * PI * t)

figure

ezplot(g(t)); xlim([-2 2]); ylim([-1.5 1.5]);

u(t) = heaviside(t)

f(t) = -t * (u(t+1) - u(t))

ezplot(f(t)); xlim([-2 2]); ylim([-1.5 1.5]);

fourier(f(t))

w(t) = u(t+(ONE/2)) - u(t-ONE/2)

fourier(w(t))

withPreferReal = simplify(ans,'Criterion','preferReal','Steps',100)

* ********************************************************************


* ********************************************************************
* Symbolic Fourier Series in Matlab (Octave)
* ********************************************************************
* Let f(t) = a0 + a1 cos(1wt) + a2 cos(2wt) + a3 cos(3wt) + ...
*               + b1 sin(1wt) + b2 sin(2wt) + b3 sin(3wt) + ...
*
*       w = 2*PI/T (where T is the period)
*
* or  f(t) = c0 + c_+1 e^(+j1wt) + c_+2 e^(+j2wt) + c_+3 e^(+j3wt) + ...
*               + c_-1 e^(-j1wt) + c_-2 e^(-j2wt) + c_-3 e^(-j3wt) + ...
*
*
* n >= 1
* a0 = (1/T) * integral(f(t), t, 0, T)
* an = (2/T) * integral(f(t)*cos(nwt), t, 0, T)
* bn = (2/T) * integral(f(t)*sin(nwt), t, 0, T)
*
* -inf < n < inf
* cn = (1/T) * integral(f(t)*e^(-jnwt), t, 0, T)
* ********************************************************************

% Define a function that is valid over one period
% For example a square wave with rising edge at t=0, period = T
% Define the function from t=0 to t=T

pkg load symbolic

syms n t T

PI = sym(pi)
ONE = sym(1)
TWO = sym(2)
ENM = sym(e)
JAY = sym(j)

u(t) = heaviside(t)



assume(T, 'positive')
assume(n, 'positive')
assumeAlso(n, 'integer')

assumptions

f(t) = u(t) - 2 * u(t - T/2)

(ONE/T)*int(f(t), 0, T) 
a0 = simplify(ans, 'Steps', 100)


(TWO/T)*int(f(t)*cos(n*(2*PI/T)*t), 0, T)
an(n) = simplify(ans, 'Steps', 100)


(TWO/T)*int(f(t)*sin(n*(2*PI/T)*t), 0, T)
bn(n) = simplify(ans, 'Steps', 100)


* ********************************************************************

(add ee301 analysis stuff)

(add plotting transform pairs)
(add plotting fourier coefficients)

(add maxima stuff for fourier series and fourier transform)

(does maxima give wrong answer for unit step?? Seems like its missina a delta term)