Contents

% /home/dimarzio/Documents/working/12229/lectures/fseries.m
% Thu Nov  2 11:27:15 2017
% Chuck DiMarzio, Northeastern University
%

Fourier series

t=[1:4096]*10e-3/4096;  % Time axis
T=1024*10e-3/4096; % Period of our function
f0=1/T
omega0=2*pi*f0;
n=[-21:21]; % range of harmonics
a=(1j./n); % Amplitudes for a square wave
a(find(n==0))=0;
a(2:2:end)=0;

f0 =
   400

Plot the coefficients

figure;stem(n,real(a));grid on;
xlabel('Freq. Index');ylabel('Real A');
figure;stem(n,imag(a));grid on;
xlabel('Freq. Index');ylabel('Imag A');

Plot them vs frequency

figure;stem(n*f0,real(a));grid on;
xlabel('f, Freq., Hz');ylabel('Real A');
figure;stem(n*f0,imag(a));grid on;
xlabel('f, Freq., Hz');ylabel('Imag A');

Sum the waves

v=zeros(size(t));
for k=1:length(n);
  v=v+a(k)*exp(1j*n(k)*omega0*t);
end;

Plot vs time. Note imaginary part is zero (it better be)

figure;plot(t,real(v));grid on;
xlabel('t, Time, sec');ylabel('real V');
figure;plot(t,imag(v));grid on;
xlabel('t, Time, sec');ylabel('imag V');

Sum the waves with a filter

n=[-101:101]; % range of harmonics
a=(1j./n); % Amplitudes for a square wave
a(find(n==0))=0;
a(2:2:end)=0;


tau=T/10;
 H=1./(1+1j*n*omega0*tau);
v=zeros(size(t));
for k=1:length(n);
  v=v+a(k).*H(k).*exp(1j*n(k)*omega0*t);
end;

Plot vs time. Note imaginary part is zero (it better be)

figure;plot(t,real(v),t,imag(v));grid on;
xlabel('t, Time, sec');ylabel('V, Volts');

What if we violate the "conjugate rule"

a=zeros(size(n));
a(98)=1;a(105)=-2j;

v=zeros(size(t));
for k=1:length(n);
  v=v+a(k).*H(k).*exp(1j*n(k)*omega0*t);
end;

figure;stem(n*f0,real(a),'b');grid on;hold on;
stem(n*f0,imag(a),'r');hold off;
xlabel('f, Freq., Hz');ylabel('A');

% Plot vs time. It's not real
figure;plot(t,real(v),t,imag(v));grid on;
xlabel('t, Time, sec');ylabel('V, Volts');

Published with MATLAB® R2017a