Fast Fourier Transform of the subtraction between a signal and the same signal but shifted
37 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Alejandro Cardona Ramirez
le 8 Mai 2019
Commenté : Greg Dionne
le 9 Mai 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')
0 commentaires
Réponse acceptée
Greg Dionne
le 8 Mai 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:
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
2 commentaires
Greg Dionne
le 9 Mai 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).
-Greg
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Transforms dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!