Fourier analysis of an spwm signal

8 vues (au cours des 30 derniers jours)
Mosharof Hossain
Mosharof Hossain le 7 Sep 2023
Commenté : Paul le 8 Sep 2023
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.

Réponse acceptée

Star Strider
Star Strider le 7 Sep 2023
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 commentaires
Star Strider
Star Strider le 7 Sep 2023
As always, my pleasure!
Paul
Paul le 8 Sep 2023
Is this line
Fs = 1/mean(diff(t));
justified given the wide variability in the delta times in the data?

Connectez-vous pour commenter.

Plus de réponses (1)

Paul
Paul le 7 Sep 2023
Modifié(e) : Paul le 8 Sep 2023
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 commentaire
Mosharof Hossain
Mosharof Hossain le 8 Sep 2023
Got it🙂 Thanks again

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by