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;
Lx = 1;
T = 1/Fs;
L = Lx*Fs;
t = (0:L-1)*T;
M = (L/2) + 1;
k = Fs*(0:(L/2))/L;
P = 2e-4 .* k.^(-3);
A = sqrt(2*P/Lx);
A(1) = 0;
A(2:end-1) = A(2:end-1)/2;
A2 = [A flip(A(2:end-1))];
phi = 2*pi*randn(1,L);
Z = A2.*exp(1j.*phi)*L;
z = ifft(Z, 'symmetric');
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;
dx = round(0.1*L)
aux1 = z(1:dx);
aux2 = z(dx+1:end);
Shiftz = [aux2 aux1];
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')
