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