Fast Fourier Transform of the subtraction between a signal and the same signal but shifted
35 views (last 30 days)
Alejandro Cardona Ramirez on 8 May 2019
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
% PSD (INPUT)
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
% SHIFTED SIGNAL
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')
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,