フィルターのクリア

Fourier analysis of an spwm signal

7 ビュー (過去 30 日間)
Mosharof Hossain
Mosharof Hossain 2023 年 9 月 7 日
コメント済み: Paul 2023 年 9 月 8 日
I have one cycle of an spwm signal as a time series data. How to do fourier analysis of the periodic signal and plot the spectrum?
The signal is as shown below.
I want to find the harmonics in this periodic signal and filter it to get back the original sine wave.
The time series spwm is attached above.

採用された回答

Star Strider
Star Strider 2023 年 9 月 7 日
Try this —
LD = load('spwm_one_cycle.mat')
LD = struct with fields:
ts1: [1×1 timeseries]
ts1 = LD.ts1
timeseries Common Properties: Name: 'unnamed' Time: [224x1 double] TimeInfo: tsdata.timemetadata Data: [224x1 double] DataInfo: tsdata.datametadata
Time = ts1.Time;
Data = ts1.Data;
figure
plot(Time, Data)
xlabel('Time (seconds)')
ylabel('Amplitude')
[FTData,Freqs] = FFT1(Data,Time);
figure
subplot(2,1,1)
plot(Freqs, abs(FTData)*2)
grid
ylabel('Magnitude')
subplot(2,1,2)
plot(Freqs, unwrap(angle(FTData)))
grid
xlabel('Frequency (Hz)')
ylabel('Phase (rad)')
sgtitle('Fourier Transform Of SPWM Signal')
function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end
.
  7 件のコメント
Star Strider
Star Strider 2023 年 9 月 7 日
As always, my pleasure!
Paul
Paul 2023 年 9 月 8 日
Is this line
Fs = 1/mean(diff(t));
justified given the wide variability in the delta times in the data?

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

その他の回答 (1 件)

Paul
Paul 2023 年 9 月 7 日
編集済み: Paul 2023 年 9 月 8 日
Hi Mosharof,
Load the data
load spwm_one_cycle.mat
Plot the differences between time samples. We see that the data is not equally spaced in time, so using discrete-time methods, like fft, is not appropriate without further massaging the data into equally spaced samples that also catch all the transitions. Maybe that could be done by using a much smaller sampling period from the original source (which looks like it might be a Simulink model?).
plot(diff(ts1.Time),'-o')
However, the signal is just the sum of rectangular pulses, and the rectangular pulse has a closed-form expression for its Continuous Time Fourier Transform (CTFT)
Here's one way we can construct the signal for analysis (this would be easier if the original source could save off the switch times and the associated signal levels at the switches).
Verify that the data has unique values:
unique(ts1.Data)
ans = 3×1
-1 0 1
Use diff to find when the signal level switches
index = [1; find(diff(ts1.Data)~=0)];
syms t w real
syms y(t) Y(w)
y(t) = 0;
Y(w) = 0;
The pattern is that the level changes only between 0 and 1, or 0 and -1. Construct the signal as the sum of rectangularPulse and build up its CTFT along the way.
for ii = 2:2:numel(index)-1
if ts1.Time(index(ii)) < 10
a = 1;
else
a = -1;
end
s(t) = a*rectangularPulse(ts1.Time(index(ii)),ts1.Time(index(ii+1)),t);
y(t) = y(t) + s(t);
Y(w) = Y(w) + fourier(s(t),t,w);
end
Verify that y(t) matches the original data
figure
subplot(211)
plot(ts1.Time,ts1.Data);
ylim([-1.1 1.1])
subplot(212)
fplot(y(t),[0 20],'MeshDensity',200)
ylim([-1.1 1.1])
fplot was taking too long to make this plot
%figure
%fplot(abs(Y(w)))
Convert to an anonymous function
Yfunc = matlabFunction(Y(w));
and evaluate over some frequency range
wval = -20:.01:20;
figure
subplot(211)
plot(wval,abs(Yfunc(wval))),grid
subplot(212)
plot(wval,180/pi*angle(Yfunc(wval))),grid
xlabel('rad/sec')
The Question stated that the data is actually one period of a periodic singal, which would be represented as a Continuous Fourier Series (CFS). The CFS coefficients are obtained by sampling the CTFT of a single period at multiples of the fundamental frequency (2*pi/T rad/sec, T = 20 sec) and scaling by 1/T. Here's a plot of the magnitude of the CFS coefficients over roughly the same frequency range used to plot the CTFT of the single period.
T = 20;
n = -60:60;
figure
stem(2*pi*n/T,(1/T)*abs(Yfunc(2*pi*n/T)))
xlabel('rad/sec')
  1 件のコメント
Mosharof Hossain
Mosharof Hossain 2023 年 9 月 8 日
Got it🙂 Thanks again

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

カテゴリ

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

Community Treasure Hunt

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

Start Hunting!

Translated by