Clear Filters
Clear Filters

Fourier Series,Spectral Analysys

7 views (last 30 days)
Tudor Valentin
Tudor Valentin on 16 May 2024
Answered: Paul on 20 May 2024
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');
  1 Comment
Paul
Paul on 16 May 2024
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.

Sign in to comment.

Answers (2)

Paul
Paul on 20 May 2024
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.

William Rose
William Rose on 16 May 2024
Edited: William Rose on 16 May 2024
[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.
  3 Comments
Tudor Valentin
Tudor Valentin on 20 May 2024
Thank you very much for the help! This looks much better than whatever i got.

Sign in to comment.

Tags

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!