Fourier Series,Spectral Analysys

Question

Tudor Valentin 2024-5-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2119626-fourier-series-spectral-analysys

回答： Paul 2024-5-20

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2024-5-20

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2119626-fourier-series-spectral-analysys#answer_1460426

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

William Rose 2024-5-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2119626-fourier-series-spectral-analysys#answer_1459091

编辑：William Rose 2024-5-16

在 MATLAB Online 中打开

@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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

William Rose 2024-5-16

在 MATLAB Online 中打开

@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 2024-5-20

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

William Rose 2024-5-20

@Tudor Valentin, you're welcome.

请先登录，再进行评论。

Fourier Series,Spectral Analysys

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（2 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Fourier Series,Spectral Analysys

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（2 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

3 个评论 显示 1更早的评论隐藏 1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论