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

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

Answer 1

[Greg Dionne]

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:

fft(z).*exp(2i*pi*dx*(0:L-1)/L)

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,

-Greg

[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.

Thanks,

Alejandro

[Greg Dionne]

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

-Greg