Fast Fourier Transform of the subtraction between a signal and the same signal but shifted

33 views (last 30 days)
My work steps are described as follows:
1. I have the Power Spectral Density PSD data which follows a power-law (in this case an equation PSD = 2e-4*k^-3, where k is frequency)
2. I convert the PSD into amplitude and assign a normal random phase to recreate a signal Z in Fourier Domain.
3. I apply IFFT to obtain the respective signal z in time domain. The size of the signal is L
4. I take this signal z and I shift it a number of elements dx. Since the signal is periodic, the "empty" space from 1 to dx is filled with the remaining signal from dx+1 to L.
5. I substract the original signal z with the new shifted signal, h.
6. I compute the FFT of this signal h.
The resultant power spectrum give me certain noise (numerical zeros) with a periodicity Dk that corresponds to L = Dk*dx.
Is there any explanation for this? Or what type of signal processing I need to perform in h?
Below, it is the code. Thanks
Fs = 1e3; %Sampling rate
Lx = 1; %Physical length of signal
T = 1/Fs; %Period
L = Lx*Fs; %Length of signal (number of points)
t = (0:L-1)*T; %Vector of time
M = (L/2) + 1; %Number of points in frequency
k = Fs*(0:(L/2))/L; %Vector of frequency
P = 2e-4 .* k.^(-3);
A = sqrt(2*P/Lx); % PSD to amplitude
% Convert to 2-side PSD
A(1) = 0; %Ensure DC offset equal zero
A(2:end-1) = A(2:end-1)/2;
A2 = [A flip(A(2:end-1))];
% Random phase shift
phi = 2*pi*randn(1,L);
% Artificial signal scaled by L
Z = A2.*exp(1j.*phi)*L;
% Inverse Fourier Transform
z = ifft(Z, 'symmetric');
% Verify Fourier Transform FFT(z) = Z
Y = fft(z);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
PSz = P1.^2*Lx./2; %Amplitude to power
dx = round(0.1*L)
aux1 = z(1:dx);
aux2 = z(dx+1:end);
Shiftz = [aux2 aux1];
%FFT of the difference between z and Shifted z
h = z - Shiftz;
Yh = fft(h);
P2h = abs(Yh/L);
P1h = P2h(1:L/2+1);
P1h(2:end-1) = 2*P1h(2:end-1);
PSh = P1h.^2*Lx./2;
loglog(k,PSz,'color','k'), hold on, loglog(k,PSh,'color','r')
xlabel('Frequency'), ylabel('Power Spectral Density')

Accepted Answer

Greg Dionne
Greg Dionne on 8 May 2019
Circularly shifting in the time domain is equivalent to multiplying by a phase ramp in the frequency domain.
The FFT of the variable "Shiftz" is equivalent to:
So the math for "Yh" is equivalent to
Yh = fft(z - Shiftz)
= fft(z) - fft(Shiftz) % FFT is a linear operator
= fft(z) - fft(z) .* exp(2i*pi*dx*(0:L-1)/L)
= fft(z) .* (1 - exp(2i*pi*dx*(0:L-1)/L))
This results in your original spectrum multiplied by a function whose zeros occur whenever exp(2i*pi*dx(0:L-1)/L) == 1.
This hopefully explains the periodicity of the zeros you are seeing.
Hope this helps,
Alejandro Cardona Ramirez
Hi Greg
This clarifies a lot the problem. Now, from the physical sense, that means those frequencies in z and ShiftZ are in phase and with the same amplitude, and therefore, the resultant value for the subtraction will vanish?
Since my plot is loglog, I will remove those points for plotting purposes.
Greg Dionne
Greg Dionne on 9 May 2019
That's right. Another way to look at it is that if z is periodic, then delaying any frequency component of z by a time interval corresponding to an integral number of periods of that frequency leaves that component essentially unchanged.
A similar effect occurs where the delay is 180 degrees out-of-phase. There the amplitude of the spectrum is doubled (instead of canceled).

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!