Piece-wise Fourier series / Inverse Fourier series

Computes the Fourier series coefficients of an arbitrarily defined piece-wise polynomial, or vice versa with visualization.

Forward Fourier series:
For function call

[R,r0,T] = fspw(c,cK,T)

Input: (see theory)
c is standard form polynomial coefficient matrix for a piece-wise polynomial
cK is the highest degree, non-piece-wise coefficient
T is the total interval measure

Output: (see theory)
R is corresponding standard form frequency domain coefficient matrix
r0 is the DC coefficient
T is the total interval measure, preserved.

Given a piece-wise polynomial function in Heaviside form defined by a real-valued coefficient matrix c_{q,k} and measured interval 0<x<T:

p(x) = \sum_{q=0}^{Q-1}[u(x-qT/Q)-u(x-(q+1)T/Q] (\sum_{k=0}^{K-1}c_{q,k}x^k) + c_{K}x^K

Function outputs the Fourier coefficients using a complex basis defined by a real-valued coefficient matrix R_{k,q}

p(x) = \sum_{n=-\infty}^{+\infty}Q_{n}e^{2\pi i n x/T}

where

Q_{n} = \begin{cases} \sum_{k=0}^{K-1}[1/(2\pi i n)]^{k+1}(\sum_{q=0}^{Q-1}R_{k,q}e^{2\pi i n q/Q}), & n \neq 0 \\ r_{0}, & n=0 \end{cases}

Inverse Fourier series:
For function call

[c,cK,T] = ifspw(R,r0,T)

Input:
R is standard form frequency domain coefficient matrix for a piece-wise polynomial
r0 is the DC coefficient
T is the total interval measure, preserved.

Output:
c is corresponding standard form polynomial coefficient matrix
cK is the highest degree, non-piece-wise coefficients
T is the total interval measure, preserved

Function output is a piece-wise polynomial function in Heaviside form defined by a real-valued coefficient matrix c_{q,k} and measured interval 0<x<T:

p(x) = c_{q,K}x^{K} + \sum_{q=0}^{Q-1}[u(x-qT/Q)-u(x-(q+1)T/Q] (\sum_{k=0}^{K-1}c_{q,k}x^{k})

Function input is the Fourier coefficients using a complex basis defined by a real-valued coefficient matrix R_{k,q}

p(x) = \sum_{n=-\infty}^{+\infty}Q_{n}e^{2\pi i n x/T}

where

Q_{n} = \begin{cases} \sum_{k=0}^{K-1}[1/(2\pi i n)]^{k+1}(\sum_{q=0}^{Q-1}R_{k,q}e^{2\pi i n q/Q}), & n \neq 0 \\ r_{0}, & n=0 \end{cases}

Plotter function:
For function call

fspwplotter(R,r0,c,cK,T)

Example 1 Forward (Stair-step Wave Fourier Expansion):
c = [1;.8;.6;.4;.2;0];
cK=0;
T = 7;
[R,r0,T] = fspw(c,cK,T);
fspwplotter(R,r0,c,cK,T)

Example 2 Inverse
R = [1 3 10 1 -4 -10; -10 4 2 5 5 5 ; -20 100 10 -10 -1 -4 ; 10 -4 -20 -30 2 3]
r0 = 0;
T = 7;
[c,cK,T] = ifspw(R,r0,T);
fspwplotter(R,r0,c,cK,T)