Fourier analysis of an spwm signal

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.

Answer 1

[Star Strider]

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

.

[Mosharof Hossain]

Hi, thank you for your answer. Can you please tell me how FFT1 calculating fourier series of the signal? More specifically I want to know how you are telling matlab that it is a periodic signal.

[Star Strider]

My pleasure!

MATLAB has no idea that it is a periodic signal, it just assumes it, as does the Fourier transform in general. The signal does not have to be periodic or infinite for the Fourier transform to work, because the Fourier transform assumes that the signal provided to it is a representative segment of an infinite (and infinitely periodic) signal.

My ‘FFT1’ function computes a one-sided fft (it assumes that the sampling periods are constant) by first calculating the sampling frequency ‘Fs’ and the associated Nyquist frequency, zero-pads the signal vector to the next power-of-2 for computational efficiency, does a windowed fft, creates a frequency vector ‘Fv’ to match the fft, and then returns them. (I created it simply because I would otherwise type all that in each time I needed to calculate a one-sided fft.) It can easily be edited to accommodate matrix signal arguments if necessary.

[Mosharof Hossain]

So fft assumes that provided signal is aperiodic (infinitely periodic) though the actual signal could be periodic. So, isn’t it going to give wrong result?

[Star Strider]

It’s not going to give a ‘wrong’ result. The result it gives is as good as the information it gets. (Using the windowing function corrects for it being finite length.) In this instance, it assumes the vector it gets as an argument is a representative section of an infinitely-repeating sequence of such segments.

[Mosharof Hossain]

Ok. Thank you very much.

[Star Strider]

As always, my pleasure!

[Paul]

Is this line

Fs = 1/mean(diff(t));

justified given the wide variability in the delta times in the data?

Answer 2

[Paul]

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')

[Mosharof Hossain]

Got it🙂 Thanks again