フィルターのクリア

Fourier Series,Spectral Analysys

19 ビュー (過去 30 日間)
Tudor Valentin
Tudor Valentin 2024 年 5 月 16 日
回答済み: Paul 2024 年 5 月 20 日 17:56
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 件のコメント
Paul
Paul 2024 年 5 月 16 日
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.

サインインしてコメントする。

回答 (2 件)

Paul
Paul 2024 年 5 月 20 日 17:56
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 2024 年 5 月 16 日
編集済み: William Rose 2024 年 5 月 16 日
[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 件のコメント
Tudor Valentin
Tudor Valentin 2024 年 5 月 20 日 11:02
Thank you very much for the help! This looks much better than whatever i got.
William Rose
William Rose 2024 年 5 月 20 日 12:07
@Tudor Valentin, you're welcome.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeSpectral Measurements についてさらに検索

タグ

製品


リリース

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by