Converting PSD to Time Domain Signal

I am currently using MATLAB and I'm thinking if I can reconstruct PSD back to time domain signal using MATLAB's ifft. I found a link in solving this problem but it's in Phyton and when I tried it in MATLAB, I encountered an error mentioning that the dimensions didn't match. This is the code that I come up with:

%PSD to time-domain signal

%periodogram using fft
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

figure(1)
plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

figure(2)
plot(t,x)
grid on
title ('Signal')

%psd to signal using ifft
magnitude = N*sqrt(psdx);
phase = 2*pi*randn(1,N);
FFT = magnitude .* exp(sqrt(-1) .* phase);
signal = ifft(FFT);

figure(3)
plot(t,signal)

I calculated the PSD using the FFT but its just so that I can do it backward from PSD to FFT and then iFFT it to find the signal. This is just an example for me to compare if doing it backwards can result in the exact same signal as the original one. But once again, MATLAB said there's an error, is there anything I can change in my code? Thank you.

• I just ran it and error is because of size mismatch. psdx, magnitude size 1x501 while phase is 1000. Also, I am curious to know if phase is random, how will you recover original signal by sqrt of magnitude alone? Mar 10 '20 at 6:50
• @jithinrj I read it in the link that its done that way, so I tried following it, I also find this link that explain about it: [link] (researchgate.net/post/…) but I'm not sure how I should code it. Mar 10 '20 at 6:54
• You can refer to this too dsp.stackexchange.com/questions/7698/… Mar 10 '20 at 7:04
• @jithinrj yes, I already link that in my questions and I already refer to it. But I chose to convert signal to PSD and then back to its signal. I tried following the Python code and change it to MATLAB version but I might miss something. Mar 10 '20 at 7:09
• As @jithinrj pointed out the data size of the phase and magnitude vector are different. This is due to xdft = xdft(1:N/2+1);, as you are taking only the half of the data points in the frequency domain. Try to change it so that you take everything, additionally you should consider to make the phase symmetrical or antisymmetrical as pointed out in the link Mar 10 '20 at 8:41

@Amelia. Hi. I tried to run your MATLAB code, but it contains several errors regarding vector lengths and produces error messages. But that's not the main issue here. The answer to your question is: If you compute the DFT of an N-point x(n) sequence you produce an N-point complex-valued xDFT(m) frequency-domain sequence. You can only reconstruct the original x(n) sequence by computing the inverse DFT of the N-point xDFT(m) sequence. You cannot recover the original N-point x(n) sequence by computing the inverse DFT of any sequence that is not the xDFT(m) sequence.

• Thank you for answering. But can we produce an xDFT(m) sequence from a PSD? Knowing that we can derive PSD from xDFT(m). I didn't know how I can do it backwards. Mar 11 '20 at 6:04
• @Amelia. Hi. No, we cannot produce an xDFT(m) sequence from a PSD. When we compute a PSD (a real-valued sequence) from an xDFT(m) (a complex-valued sequence) the PSD does not contain the phase information of the xDFT(m) sequence. For example, if I have the complex number 1 +j2 I can compute its magnitude to be roughly equal to 2.2361. Now if I e-mailed you and wrote, "Amelia, I have a complex number whose magnitude is 2.2361. Please tell me what are the real and imaginary parts of my complex number." You would not be able to tell me what are the real and imaginary parts of my complex number. Mar 11 '20 at 10:47

At the risk of sounding stupid, I am going to give a fresh perspective (different from the link referenced above) to this problem. I will list my reasoning and finally give the code I tried.

1. When the signal length is 1000, auto-correlation will have 1999 length. Hence FFT of PSD should have at least 1999 length. In my example I used 2000 for simplicity.
2. I did not understand why multiply by 2 is needed, so I am going to skip that.
3. PSD is symmetric (since auto-correlation is real valued) but only half the frequencies are used. I feel this is incorrect. While plotting this is fine. But while taking IFFT, what we need to do is to take 1:2:2000 instead of taking first 1000 values of PSD.
4. Although the input signal is real, due to random phase addition, the final signal is complex. So I have plotted only the real part.

%PSD to time-domain signal
%periodogram using fft
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));

N = length(x); xdft = fft(x,2000); %xdft = xdft(1:N/2+1); psdx = (1/(Fs*N)) * abs(xdft).^2; %psdx(2:end-1) = 2*psdx(2:end-1); freq = 0:Fs/length(xdft):Fs*(1-1/length(xdft));

figure(1) plot(freq,10*log10(psdx)) grid on title('Periodogram Using FFT') xlabel('Frequency (Hz)') ylabel('Power/Frequency (dB/Hz)')

figure(2) plot(t,x) grid on title ('Signal')

%psd to signal using ifft magnitude = N*sqrt(psdx(1:2:end)); phase = 2*pi * rand(1,N); FFT = magnitude .* exp(sqrt(-1) .* phase); signal = ifft(FFT);

figure(3) plot(real(signal))

• I tried your code and I realized that the ifft of FFT is still in a complex numbers that's why the final signal is complex. But why is it when I ifft the fft(x) it will not result to a complex variable? Is there something wrong with the FFT equation? Mar 11 '20 at 6:44
• @AmeliaLita you will get a complex signal. Because, your original signal (|x|.^2) was real, the magnitude is symmetric and phase is anti-symmetric. By adding random phase, it is no longer anti-symmetric. Hence the resulting IFFT sequence will no longer be real. Try making the phase anit-symmetric random values. You should get real signal. Mar 11 '20 at 10:46
• How can I make an anti-symmetric random values for the phase in MATLAB? I'm sorry if this sounds stupid. Mar 12 '20 at 16:01
• @Amelia. Hi. If you still need help then please e-mail me at: R.Lyons@ieee.org. I will help you. Mar 13 '20 at 9:46

Try the following code. Note that you need to make the second half of the phase data negative of the first half put in reverse order, such that the second FFT half is a complex conjugate mirror of the first half.

clear
clc

% time history
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = randn(size(t));

% one-sided periodogram using fft, rectangular window
N = length(x); % N should be an even integer to make the Nyquist frequency as one of the FFT entries
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

plot(freq,10*log10(psdx), 'LineWidth', 1, 'DisplayName','PSD from Original Time-History')
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

% one-sided psd to signal using ifft and Fs
Fs=Fs;
psdx(2:end-1) = 0.5 * psdx(2:end-1);
psdx = [psdx flip(psdx(2:end-1))]; % converting one-sided psd to two-sided
fft_mag = sqrt(N*Fs*psdx);
fft_phase = 2*pi*randn(1,N/2+1);
fft_phase = [fft_phase -flip(fft_phase(2:end-1))];
FFT = fft_mag .* exp(sqrt(-1) .* fft_phase);
signal = ifft(FFT);
signal = real(signal); % you can check that phase angles are very small

% one-sided periodogram using fft, rectangular window
N = length(signal);
xdft = fft(signal);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(signal):Fs/2;

hold on
plot(freq,10*log10(psdx), '--r', 'LineWidth', 2, 'DisplayName','PSD from Converted Time-History')
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
legend
• Welcome to SE.SP! Please just use the {} button to format code. No need to insert HTML tags to get it formatted.
– Peter K.
Nov 4 '21 at 20:33