Does anyone know how to construct original discrete time signal by using IFFT(x) function?

Hi,

I need to construct my 24 point discrete time series by using its first 5 harmonics. I know I can use x=ifft(y) to obtain the original signal in which x and y should have the same dimensions (24x1 and 24x1). My question is that is it possible to obtain the same size of x (24x1) by using let's say first 5 element of y?

Thanks,

[Adam]

24 samples doesn't seem very many to expect to have 5 distinct harmonics. But if you know their frequencies and they fit on the grid defined by just 24 samples in frequency space, then yes you can recreate your signal if you have the sampling frequency.

Answer 1

[Star Strider]

You could use only the first five frequencies, if you wanted an approximation. Reconstructing from less than all the available frequency components loses detail in the reconstruction. The reconstructed signal will be the size of the frequencies you used to reconstruct it. You can zero-pad the reconstructed signal out to the length of the original vector (I did not do that here), but you cannot recover the information you lost.

Example:

t = linspace(0, 1, 24);                                     % Time Vector
x = (t >= 1/3) & (t <= 2/3);                                % Signal Vector

y = fft(x);                                                 % Fourier Transform
ys1 = fftshift(y);                                          % Shift To Centre
ys2 = fliplr(fftshift(ys1(8:18)));                          % Take Centre 5 Frequencies, Shift Back, And Flip

xinv = ifft(ys2);

figure(1)
plot(x,'b')
hold on
plot(xinv,'r')
plot(ifft(y), '--g')
hold off
legend('Original Waveform', 'Reconstructed From First 5 Freqencies', 'Reconstructed From All Frequencies')

[Fred]

Thanks for the answer. I have another question about the use of zero-pad in frequency domain. When i tried to use

xinv = ifft([0 0 0 0 0 0 0 ys2 0 0 0 0 0 0]);

instead of using

xinv = ifft(ys2);

I got so weird reconstructed signal. Can you just help me how to use zero-padding in the frequency domain so that I can get my original signal with a small error margin?

Thanks

[Star Strider]

My pleasure.

The fftshift function in ‘ys1’ rotates ‘y’ by half the length so that 0 Hz is in the centre of the vector rather than at the beginning. In ‘ys2’ with your 24-element vector (and the 24 element fft it creates), the index range of 8:18 of the shifted vector incorporates the 5 frequencies to the ‘left’ of centre, the centre 0 Hz, and the 5 frequencies to the ‘right’ of centre, since symmetry is important. I then shifted it back, and discovering that it didn’t shift the way I wanted it to (so 0 Hz is the first element) use fliplr to get it in the correct orientation. Then I did the ifft on it to get the result.

My code wouldn’t work if you don’t use those two assignments!

The demonstration would have been more impressive with a longer data vector, because I then would have been able to demonstrate the Gibbs phenomenon. It’s definitely visible, but difficult to explain in the plot with so few data. It’s responsible for the ‘zig-zag’ appearance in the red line, due to the loss of high frequency components.

-------------------------------------------------------

EDIT — Don’t zero-pad in the frequency domain. The purpose of zero-padding a time-domain signal is to increase the frequency resolution in the transformed spectrum. A time-domain signal is not necessarily symmetrical (there’s no reason to expect it to be). The frequency domain signal is symmetrical, making zero-padding, even with symmetrical zero vectors, problematic. I’ve never zero-padded a frequency domain signal, and I can’t find any discussion of it in my signal processing references (quick search). If you can find such a reference, please share it. I experimented with this just now, and I can’t recover anything useful.

It apparently it is possible to interpolate using the fft, and the interpft function does just that. The documentation doesn’t go into any detail as to how it works.

[Fred]

Thanks for the reply.

How can i reconstruct my whole 24 points discrete x(n) signal without using zero-padding then?

I have a paper in which it says you can reconstruct your daily load profile (24 discrete points each represents average load value at corresponding hour, e.g, 1, 2, 3,...,24 hr) by using only first 4 harmonic components.

If there is no way to get 24-point discrete signal by using ifft(y1), where y1=y(1:4) (first 4 frequency components of the 24 component vector of y, where y=fft(x)), how these guys are able to get original signal? Do you have any idea about this?

[Star Strider]

My pleasure.

You would have to attach a PDF of the paper here (use the ‘paperclip’ icon) so that I could read it to understand what they are doing and their methods. I probably don’t have free access to it.

They might simply be using a discrete lowpass filter to pass only the first 4 frequencies in the time-domain signal. That’s how I would do it. (The Nyquist frequency would be 1/12 cycles/hour, so the normalised frequency for the filter would be 4/12 or (1/3)*pi radians.) A FIR filter or a Type II Chebyshev design would be my choice. It will be interesting to see how those authors do it.

[Fred]

Here is the paper that I'm talking about.

[Star Strider]

That’s an interesting concept.

If you want to do the frequency domain analysis, with only 24 points, you can easily use a loop to calculate Eqns (4) and (6). I doubt that there would be anything to be gained by using the fft function, considering the data you want to recover and the length of the records only being 24 samples.

That said, this seems to be taking the long way round to avoid using a simple filter. The only advantage I can see to this approach is if the filters to isolate the various harmonics are longer than the data record. That would be the situation with a FIR filter, but might not be with an IIR filter. The paper didn’t discuss the advantages of their approach over using a discrete filter, and since I have no experience in doing load forecasting, I’ve not dealt with the practical aspects of designing filters for such data.

See if this filtering approach does what you want:

The Code —

t = linspace(0, 1, 24);                                     % Time Vector
x = (t >= 1/3) & (t <= 2/3);                                % Signal Vector

Fs = 24;                                            % Sampling Frequency (Samples/Day)
Fn = Fs/2;                                          % Nyquist Frequency
Wp = [0.6  1.4]/Fn;                                 % First Harmonic Passband
Ws = [0.2  1.8]/Fn;                                 % First Harmonic Stopband
Rp = 10;
Rs = 30;
[n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);
[z,p,k] = cheby2(n,Rs,Ws);
[sos1,g1] = zp2sos(z,p,k);                          % Filter For First Harmonic

figure(1)
freqz(sos1, 2^16, Fs)

Wp = (1+[0.6  1.4])/Fn;                             % Second Harmonic Passband
Ws = (1+[0.2  1.8])/Fn;                             % Second Harmonic Stopband
[n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);
[z,p,k] = cheby2(n,Rs,Ws);
[sos2,g2] = zp2sos(z,p,k);                          % Filter For Second Harmonic

figure(2)
freqz(sos2, 2^16, Fs)

H1 = filtfilt(sos1,g1,+x);
H2 = filtfilt(sos2,g2,+x);

figure(3)
plot(t,x, 'LineWidth',1)
hold on
plot(t, H1)
plot(t, H2)
hold off
grid
set(gca, 'YLim',[-0.5  1.5])
legend('Original Signal', '1^{st} Harmonic', '2^{nd} Harmonic')

The Bode plots are stated in units of Hz. In reality, they are in units of ‘cycles/day’.

The ‘+’ in front of the ‘x’ in ‘+x’ converts the logical vector to a numeric vector. It would not make any difference with real data.

The Plot —

[Fred]

Thanks Star, I'll check it out and let you know.

[Star Strider]

My pleasure.

Note that all I did to design the second harmonic filter was to add 1 to the bandpass and bandstop vectors. Do the same for as many harmonics (up to 12) as you want. The Chebyshev Type II design will provide good separation with filters short enough to work with a 24-element data vector.

Answer 2

[Lugo Hino Far]

Help! I need to build those charts. They are Fourier Series.