Fourier Series,Spectral Analysys

I am trying to implement the spectral analysys for the signal represented in Fig.8. Im trying to do this by using the Fourier Series to generate the signal(it doesnt have to be perfect like in Fig.8) however I am very lost at this point and I dont think thats how the spectral analysys should be looking,i kept trying but I am not able to do it.

Thank you very much for your time.

Sorry for the words not being in english

% Parametrii semnalului dreptunghiular
t1 = 0.2; % Momentul de pornire al creșterii
t2 = 0.8; % Momentul de sfârșit al creșterii
A = 1; % Amplitudinea semnalului dreptunghiular

% Perioada semnalului dreptunghiular (este 1 secundă în acest caz)
T = 1;

% Numărul de armonici pentru seria Fourier (poți ajusta acest număr pentru o aproximare mai bună)
N = 10; % Vom folosi 10 termeni (5 perechi de cosinus și sinus)

% Vectorul de timp
t = linspace(0, T, 1000); % 1000 de puncte între 0 și T pentru o reprezentare mai precisă

% Inițializare semnal dreptunghiular aproximativ
f = zeros(size(t));

% Calcul coeficienți Fourier
a0 = 1; % Coeficientul mediei (pentru semnalul dreptunghiular cu amplitudinea 1)

for n = 1:N
    freq = 2*n - 1; % Frecvența armonicului (doar frecvențele impare)
    an = (2*A/T) * (1/(pi*freq)) * (sin(2*pi*freq*t2/T) - sin(2*pi*freq*t1/T)); % Coeficientul an
    bn = 0; % Pentru semnalul dreptunghiular, coeficientul bn este 0 (nu există termeni sinusoidal)

    % Adăugare termen în sumă
    f = f + an * cos(2*pi*freq*t);
end

% Adăugare termenul de medie
f = f + a0/2;

% Plotare semnal dreptunghiular aproximativ
figure;
plot(t, f, 'b', 'LineWidth', 2);
xlabel('Timp [secunde]');
ylabel('Amplitudine');
title('Semnal Dreptunghiular Aproximat prin Seria Fourier');
grid on;

%Analiza spectrala
Y=fft(f);
Pyy=Y.*conj(Y);
s=length(Pyy);
f=20000*(0:s/2-1)/s;
figure;
plot(f,abs(Pyy(1:s/2)));
xlabel('Frecventa [Hz]');
ylabel('Spectru'); 

[Paul]

Hi Tudor,

I don't think that the expressions for a(n), b(n), or a0 are correct. If those expressions were developed using Matlab, consider posting that code as well.

Answer 1

[Paul]

Define the function over one period

T0 = 1; % period
f0 = 1/T0;
syms t
s(t) = rectangularPulse(.2,0.8,t);
figure
fplot(s(t),[0 1]);axis padded
xlabel('t');ylabel('s')

The sine-cosine Fourier series coefficients are:

a = @(n) -(sin((2*pi*n)/5) - sin((8*pi*n)/5))./(n*pi);  % n >= 1
b = @(n)  0*n; % n >= 1
a0 = 6/5;

Reconstruct the signal using the first 15 harmonics and add to plot

t = 0:.01:T0;
xv = 0*t + a0/2;
for k = 1:15
    xv = xv + a(k)*cos(k*2*pi*f0*t) + b(k)*sin(k*2*pi*f0*t);
end

hold on
plot(t,xv)

To approximate the continuous Fourier series coefficients with the discrete Fourier series coefficients as computed by fft we need to first sample the signal. We need to enusre that the sampling period, Ts, is selected such that T0/Ts is an integer and that integer must be at least 31 if we want to compare to the first 15 harmonics.

Ts = 1/31;
N  = T0/Ts

N = 31

tval = (0:N-1)*Ts;
sval = double(s(tval));
stem(tval,sval)

Note that this selection of sampling period doesn't sample on the edges of the pulse.

Compute the discrete Fourier series coefficients and scale by Ts to compare to continuous Fourier series coefficients.

c = fft(sval)*Ts;

Compare at f = 0

[c(1) a0/2]

ans = 1x2

0.5806 0.6000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Because s(t) is real, we can take advantage of the symmetry in c to compute a and b as follows.

figure
stem(1:15,a(1:15)),hold on;stem(1:15,2*real(c(2:16)))
xlabel('n');ylabel('a(n)')

figure
stem(1:15,b(1:15)),hold on,stem(1:15,2*imag(c(2:16)))
xlabel('n');ylabel('b(n)')

Play around with higher values of Ts to see how the results change.

Answer 2

[William Rose]

@Tudor Valentin,

[edit: fix spelling errors]

Here is an example to assist you. The example below uses fft() to compute the spectrum. It displays the full spectrum and the first N harmonics. You can compute an and bn using the standard formulas. Your results with an and bn should look similar to the results below. There may be a difference of a constant scaling factor, depending on how the coefficients an and bn are normalized.

Define constants:

A=1; T=1; t1=0.2; t2=0.8;
fs=1000;            % sampling frequency (Hz)
N=10;               % number of harmonics

Create the time vector and the signal:

N1=T*fs;            % number of time values
t=(0:(N1-1))/fs;    % time vector
x=zeros(1,N1);
x(t>=t1 & t<t2)=A;  % x(t)
% Plot x versus time
plot(t,x,'-r.');
xlabel('Time (s)'); ylabel('Amplitude'); title('x(t)')

Find power spectrum using fft() function. Subtract the mean value from x first. This is a standard step when computing the power spectrum.

X=fft(x-mean(x));       % x and X are different
Pxx=abs(X(1:N1/2+1));   % power spectrum
f=fs*(0:N1/2)/N1;       % vector of frequencies

Plot power spectrum. Second plot is same as first plot, but second plot shows only the first N harmonics.

figure;
subplot(211), plot(f,Pxx,'-b.');
grid on; xlabel('Frequency (Hz)'); ylabel('|X(f)|^2'); title('Power Spectrum')
subplot(212), plot(f(1:N+1),Pxx(1:N+1),'-b.');
grid on; xlabel('Frequency (Hz)'); ylabel('|X(f)|^2'); title('Power Spectrum')

Good luck with your work.

[William Rose]

@Tudor Valentin,

I now understand that you wish to perform a Fourier synthesis using the first N harmonics.

You used a formula that has bn=0 for all n. Therefore I will do the same. This will produce an even function, since each term in the sum is a cosine function. Therefore the resulting function x(t) will be shifted in time relative to the funciton in your original post. It will have the correct pulse width and correct amplitude and mean value.

where

a0=A*(t2-t1)/T = mean value of the signal

and

Define constants.

A=1;                        % amplitude
T=1; t1=0.2; t2=0.8;        % times (s)
Nh=10;                      % number of harmonics
N=1000;                     % number of points in the signal

Define time vector.

dt=T/N;                     % time step size (s)
t=(0:N-1)*dt;               % time vector (s)

Compute x(t).

w0=2*pi/T;                  % fundamental frequency (rad/s)
a0=A*(t2-t1)/T;
x=a0*ones(1,N);
for n=1:Nh
    a=(2*A/(n*pi))*sin(n*pi*(t2-t1)/T);
    x=x+a*cos(n*w0*t);
end

Plot x(t).

plot(t,x,'-r.');
xlabel('Time (s)'); ylabel('Amplitude'); grid on

The function x(t) above is ~=1 for 60% of the cycle, and ~=0 for 40%. This matches your original signal.

[Tudor Valentin]

Thank you very much for the help! This looks much better than whatever i got.

[William Rose]

@Tudor Valentin, you're welcome.