Signals&Systems Laboratory 2013/2014
Experiment 3.
Fourier Series. Signal expansion and synthesis.
Introduction

Fourier series.
Consider a periodic signal x(t) with period T. If some conditions (see Dirichlet conditions from lecture book) are satisfied, the function (signal) x(t) can be expressed as series:
(1)
where .
Coefficients of the trigonometric Fourier series may be calculated from the following relations:
(2)
(3)
(4)
where t_{0} is any real number.
Formula (1) can be also expressed in more convenient equivalent form:
(5)
Cosine terms in (1) and (5) are called harmonic (nth harmonic has amplitude, phase and angular frequency , constant term (2) is the 0th harmonic. First harmonic:
(6)
is called usually basic harmonic with angular (basic) frequency . The plot of amplitude A_{n} i terms of n is called amplitude spectrum and plot of phase as a function of n is called phase spectrum. Amplitudes and phases may be calculated from the relations:
(7)
(8)
The periodic signals are often presented in more compact form called exponential Fourier series:
(9)
where
. (10)
Relations between different Fourier series forms are specified as follows:
. (11)
Similarly, the plot of versus n is called the amplitude spectrum and the plot of versus n is called phase spectrum. Note, that following relationships are valid:

How to compute coefficients.
The amplitude and phase spectra, based on formulas (5) or (10) can be computed by use of Matlab integration formulas like:
int
Symbolic integration
Syntax
int(expr,var,a,b) computes the definite integral of expr with respect to var from a to b.
More often, Fourier series computations are performed numerically using the Discrete Fourier Transform (DFT), which in turn is implemented numerically using an efficient algorithm known as the Fast Fourier Transform (FFT). Assuming, that sampling frequency used for discretization of the periodic function interval is equal to f_{s} (samples per second) and applied DFT transform has Nelements, basic harmonic of the Fourier series may be calculated from he formula:
. (12)
Frequence of the next stripes of Fourier series spectrum are of the form:
. (13)
Constant term and amplitudes of the harmonics (up to i=N/2) of the Fourier series having form(5) follows the relations:
(14)
Phases of the harmonic are computed from the arguments (command angle)of X_{n}._{ }
_{ }
As an illustration, the following code (core of the mfile) shows how to use fft approach to obtain Fourier expansion coefficients. You can study this code and further enhance it to complete your work.
function [SpecAmp,SpecPh,Y1]=Spectrum_by_FFT_2013(fun,N,T)
%N==> dlugosc transformaty DFT (Length of DFT)
%T==> czas pomiaru (wielokrotnosc okresu?) (oservation time, suggested:multiplicity of the period T)
ts=(T/N);
fs=1/ts;%
t1=0:ts:T;
ta=0:0.01*ts:T;
n1=0:1:N/2;
y1=feval(fun,t1); %function 'fun' after disretization for DFT
ya=feval(fun,ta); %function 'fun' after discretization enabling to print analog curve
Y1=fft(y1,N); %
Y1(2+N/2:N)=[]; %
SpecAmp=abs(Y1)*2/N; %
SpecAmp(1,1)=SpecAmp(1,1)/2; %
SpecPh=(abs(Y1)>0.01).*angle(Y1)*180/pi;
y=synth(T,N,SpecAmp,SpecPh,ta);
subplot(5,1,1);plot(ta,ya);grid on;
subplot(5,1,2);stem(n1,SpecAmp);grid on;
subplot(5,1,3);stem(n1,SpecPh);grid on;
subplot(5,1,4:5);plot(ta,y);grid on;
function y=synth(T,N,Amp,faza,x); %
y=zeros(1,length(x));
for i=1:N/2
y=y+Amp(i+1).*(cos(i*2*pi*(T^1).*x+faza(i+1)*pi/180));
end
y=y+Amp(1);
return
Useful Matlab functions: fft(x,N), length(), abs, angle, stem figure, xlabel, ylabel, title.
Numerical calculations may be performed also by Simulink models (available only in Matlab 5.3 installation). Exemplary scheme is shown in Fig.1
Fig.1. FFT based algorithm for finding harmonic’s coefficients – Simulink implemetation.
Fig.3.Block parameters..

Laboratory work
For given set of periodic signals calculate Fourier series coefficients (for all three forms) using

FFTbased method  Matlab script (open Matlab 5.3) delivered by your teacher (from website).

Numerical integration (Symbolic Toolbox from Matlab 2013b)

Find Fourier expansions (approximations) in the form (1) , (5) and (9). Complete Table with coefficients, plot spectra. Discuss the coefficients of various Fourier series in relation to the symmetry of the signal.

Observe the Gibbs phenomena (changing number of synthesized harmonics). Briefly explain!

Observe, what happens, if the sampling period T_{s} is not an integer multiplicity of the signal period T.
Fig.4a
Fig.4b
Fig.5
Exemplary mfile enabling us to calculate (and plot) Canonical Fourier Series Coefficients
% Plotting the Fourier Series Coefficients
% Współczynniki kanonicznej postaci Szeregu Fouriera
syms t k n
%x=exp(t);
x=triangle1(t); % funkcja zawierająca opis sygnału (function containing signal description)
to=0;
T=4;
w=2*pi/T;
ao=(1/T)*int(x,t,to,to+T);
a=(2/T)*int(x*cos(n*w*t),t,to,to+T);
n1=1:6;
a=subs(a,n,n1);
an=[eval(ao) a];
n2=[0 n1];
b=(2/T)*int(x*sin(n*w*t),t,to,to+T);
n1=0:6;
bn=subs(b,n,n1);
subplot(2,1,1);
A=sqrt(an.^2+bn.^2);
stem(n2,A);grid on
title('Canonical Fourier Series Coefficients: Amplitude An')
legend('An,n=0:6')
subplot(2,1,2);
thita2=atan2(bn.*(abs(bn)>0.0001),an)*180/(1*pi);
stem(n1,thita2);grid on
title('Canonical Fourier Series Coefficients: Phase theta_n')
legend('\angle thita,k=0:6')
