0
$\begingroup$

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.

| improve this question | | | | |
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – jithin Mar 10 at 6:50
  • $\begingroup$ @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. $\endgroup$ – Amelia Lita Mar 10 at 6:54
  • $\begingroup$ You can refer to this too dsp.stackexchange.com/questions/7698/… $\endgroup$ – jithin Mar 10 at 7:04
  • $\begingroup$ @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. $\endgroup$ – Amelia Lita Mar 10 at 7:09
  • $\begingroup$ 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 $\endgroup$ – Irreducible Mar 10 at 8:41
1
$\begingroup$

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

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Amelia Lita Mar 11 at 6:04
  • $\begingroup$ @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. $\endgroup$ – Richard Lyons Mar 11 at 10:47
0
$\begingroup$

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

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Amelia Lita Mar 11 at 6:44
  • $\begingroup$ @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. $\endgroup$ – jithin Mar 11 at 10:46
  • $\begingroup$ How can I make an anti-symmetric random values for the phase in MATLAB? I'm sorry if this sounds stupid. $\endgroup$ – Amelia Lita Mar 12 at 16:01
  • $\begingroup$ @Amelia. Hi. If you still need help then please e-mail me at: R.Lyons@ieee.org. I will help you. $\endgroup$ – Richard Lyons Mar 13 at 9:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.