How to do a phase shift of a signal from FFT output data

Hello,

I'm dealing with the post-processing of the simulation results of a BLDC motor Simulink model operating with a PWM control method. My goal was to export the phase current and analyze it in Matlab to get the FFT parameters (magnitudes of the different harmonics) in order to use the signal afterward by means of these coefficients (according to the definition of the FFT coefficients) -in other words, I need to use the current in another software to estimate additional losses due to current harmonics-. I was able to import a frame of the current waveform corresponding to a fundamental period and doing the FFT:(here is the original signal imported from Simulink)

and here the signal built by means of the FFT coefficients, with its fundamental:

the aim now would be to shift the whole signal in order to let it be in phase with respect to cos(x); here is the piece of code I've been writing so far:

N = 14970;
x = Current_per.Data; % vector of current from Simulink result
h = ((1:N/2) - 1)'; % harmonic orders
[ah,bh,X]=fft_M(x);
x1=linspace(0,2*pi,N/2);
I1 = ah'*cos(h*x1) +bh'*sin(h*x1); % current reconstruction with fft data
figure
hold on
plot(x1,I1)
plot(x1,ah(2)*cos(x1) + bh(2)*sin(x1))
% plot(x1,ah(6)*cos(5*x1) + bh(6)*sin(5*x1))
figure
bar(h,sqrt(ah.^2+bh.^2))
axis([-10 500 -0.01 0.05])

I attach the .mat file of the Current_per data; and here is the implementation of the fft_M:

function [ah,bh,X]=fft_M(x)

% remember to drop the last point of x2 if equal
% to the first one before passing it to fft

N = length(x); % number of samples
if x(N)==x(1)
    x=x(1:N-1);
end
X = 2/N*fft(x); % make the complex FFT
X = X(1:floor(N/2)); % drop the last half
X(1) = X(1)/2; % mean value correction
X(2:end) = X(2:end).*exp(-1i*angle(X(2)));
ah = real(X); % cos harmonic coefficients
bh = -imag(X); % sin harmonic coefficients
end

My idea, as you can notice, was to shift the all harmonics by the angle of the fundamental in order to simply move the above wave, but the result is not only a simple shifting:

I do really hope I am blind to a possible mistake and that someone could help me; please notice that I have already considered the possibility to take directly a current wave which is in phase with cos(x) but it's not easy and furthermore, I may need to shift the signal afterward anyway (so I need to fix the fact that with the shifting method I use, the signal is no longer the original one). Any help would be really appreciated. Thanks in advance

Answer 1

[David Goodmanson]

Hi matteo,

MODIFIED

here is an approximate way to do this. The code below stays with the entire frequency array, rather than cutting it in half, taking real & imaginary parts, multiplying by 2 etc, which I find to be prone to error, at least for me. I did use fftshift to shift zero frequency to the center of the corresponding frequency array. N is assumed to be even.

Shifting the time domain signal x(t) to x(t-t0) involves multiplying the fourier transform by a factor

exp(-2*pi*i*f*t0)         [1]

in the frequency domain.

If y = fft(x), then the contribution of the fundamental is

y(ip)*exp(2*pi*i*f(ip)*t) + y(im)*exp(2*pi*i*f(im)*t)

where y(ip) is the amplitude at the the first fundamental positive frequency f(ip), and y(im) is the amplitude at the first fundamental negative frequency f(im) = -f(ip). For a real signal x(t), y(ip) and y(im) are complex conjugates. To get a pure cosine wave with no sine contribution, y(ip) and y(im) have to be real. So one can find the phase of y(ip) and multiply the y array by [1] in such a way that the phase factor

exp(-2*pi*i*f(ip)*t0)

eliminates the phase in y(ip) and also y(im). Then transform back. To obtain a negative cosine wave which is close to what you started with, you can add pi to angle(y(ip)).

The fft is associated with a time array with equally spaced points, with spacing denoted by delt. The problem with the method is that things work out perfectly only if t0 is an exact multiple of delt (so that you are translating the waveform by an exact number of time steps). That is hardly ever going to be the case, and if it is not then you get an approximate solution that is probably not too bad if N is large. But even then it can have problems.

In the code I arbitrarily assumed a sampling frequency of 1 kHz although the value chosen for fs does not affect the result.

x = Current_per.Data';              % transpose
n = 14970;
fs = 1000;
delt = 1/fs;
delf = fs/n;
t = (0:n-1)*delt;
f = (-n/2:n/2-1)*delf;
y = fftshift(fft(x))/n;
i0 = n/2+1;
ip = i0+1;
im = i0-1;
first_fund = y(ip)*exp(2*pi*i*f(ip)*t) + y(im)*exp(2*pi*i*f(im)*t);
figure(4)
plot(t,x,t,first_fund)              % same as what you have

t0 = (angle(y(ip))/(2*pi*f(ip)))
y1 = y.*exp(-2*pi*i*f*t0);          % y1(ip) and y1(im) are real.
y1(1) = real(y1(1));                % band aid fix                
x1 = n*(ifft(ifftshift(y1)));       % approx. circularly shifted signal
first_fund1 = y1(ip)*exp(2*pi*i*f(ip)*t) + y1(im)*exp(2*pi*i*f(im)*t);
figure (5)
plot(t,x1,t,first_fund1)

In the shifted y array, y(1) is the nyquist frequency. Aside from the dc term it is the only frequency that does not pair up with another frequency in the [positive f, negative f] sense. To obtain a real result when going back to the time domain, this term should be real, but when t0 is not an exact multiple of delt it is not. The band aid fix is an approximation that removes the imaginary part.

[matteo leandro]

Thorough and effective answer; problem solved!! Thank you so much!

[David Goodmanson]

You're welcome. One thing I should have been more clear on is the use of 'real' to eliminate the small imaginary part. That meant eliminating small nuisance terms down around 10^-8i (I was surprised they were not even smaller), as opposed to ignoring .01i or something like that which would have meant that the code wasn't working.

[matteo leandro]

Good to know it! Thank you again then! I'm now taking advantage of your willingness, but I have another question: is there a way to get the value of the signal pointwise in the time domain by using the transformed signal, or something similar? I am currently using the FFT coefficients for the positive frequencies to rebuild the time domain signal (with the same coefficients I am able to evaluate the signal value pointwise); however, it does not seem to be accurate since the signal rebuilt in this fashion looks a bit different from the one I get by simply anti-transforming it (I have also noticed that it goes better if I increase the number of samples, but it also takes more time -of course-).

[David Goodmanson]

Hi matteo,

I modified the answer after realizing that the method is inherently less exact than I thought. The 10^-8 issue in my previous comment goes away but only by acknowledging the approximate nature of the method.

As far as getting the signal back pointwise in the time domain, I thought that is what the code does, although it's only going to be exact when t0 is an exact multiple of delt (see modified answer). One possibility is to choose t0 to be the closest exact multiple of delt. Then the ffts work out perfectly although then the fundamental cos term has a small phase shift.

One thing worth emphasizing is that the translation expression

exp(-2*pi*i*f*t0)         [1]

depends on f, and multiplies every term in the frequency data array by a different phase factor. Your original approach multiplies every term in the array by the same phase factor which is not correct.