1

I would like to generate some figures to show the disadvantages of using Fourier transforms for time series analysis. I aim to show that 2 apparently very distinct signals have a spectra with very similar shape. To begin I create my series:

t =  0:.01:2;
y = sin(2.*pi.*5.*t)+sin(2.*pi.*10.*t);
r1 = find(t <= 1);
r2 = find(t > 1);
y2a =  sin(2.*pi.*5.*t(r1));
y2b = sin(2.*pi.*10.*t(r2));
y2 = [y2a,y2b];

figure(1);
subplot(211);
plot(t,y,'k');
subplot(212);
plot(t,y2,'k');

Producing: enter image description here

Now, I would like to show that their spectra have very similar shapes: enter image description here

These are examples taken from some class notes that I would like to reproduce in matlab. However, I am having difficulties in reproducing the second plot. Could anyone suggest how the second plot can be produced in matlab with the information provided?

Eitan T
  • 32,660
  • 14
  • 72
  • 109
KatyB
  • 3,920
  • 7
  • 42
  • 72
  • have you tried `doc fft` ? this will also give you more than a hint: http://stackoverflow.com/questions/10758315/understanding-matlab-fft-example – bla Dec 09 '12 at 09:26
  • @Kate this is not the issue, but `t` should be from 0 to 4, according to your class notes, but you have defined it from 0 to 2. – Eitan T Dec 09 '12 at 09:33
  • Ok, so saying that I change it from 0 to 4, how would I reproduce the second plot? – KatyB Dec 09 '12 at 09:48

1 Answers1

2

Reproducing those plots is fairly easy. However, note several things:

  1. In your class notes, the time domain plots are from t = 0 to t = 4 (not t = 2 like you did).
  2. The frequency domain plots show only the positive frequencies, i.e. half the frequency spectrum. The frequency range is from 0 to 20Hz, meaning that the sampling frequency is 40Hz.
  3. The amount of energy must remain the same, so if only half the spectrum is displayed, the plotted frequency components should be multiplied by a factor of 2.

That said, here's the complete code (including fixes to your time domain plots):

% # Time domain plots
fs = 40;
t = 0:(1 / fs):4;
y1 = sin(2 * pi * 5 * t)+ sin(2 * pi * 10 * t);
y2 = [sin(2 * pi * 5 * t(t <= 2)), sin(2 * pi * 10 * t(t > 2))];

figure
subplot(2, 1, 1), plot(t, y1)
subplot(2, 1, 2), plot(t, y2)

% # Frequency domain plots
Fy1 = abs(ifft(y1));
Fy2 = abs(ifft(y2));
N = numel(t);
idx = 1:numel(Fy1) / 2;    % # Indices of half the spectrum
f = fs * (0:(N - 1)) / N;  % # Actual frequencies

figure
subplot(2, 1, 1), plot(f(idx), 2 * Fy1(idx))
subplot(2, 1, 2), plot(f(idx), 2 * Fy2(idx))

The time domain plots are:

enter image description here

and the corresponding frequency domain plots are:

enter image description here

Eitan T
  • 32,660
  • 14
  • 72
  • 109
  • 1
    Great, thank you very much. I used the periodogram function to compute the power spectra but like you said the plots only show the positive frequencies, this caused some confusion. – KatyB Dec 09 '12 at 10:49
  • @Kate Ah, I have found a mistake in my frequency plot: the plot shows only half the spectrum, hence everything (including the peaks) must be twice as high, in order to maintain the same amount of energy. Therefore, I've amended my answer by multiplying the plotted `Fy1` and `Fy2`. – Eitan T Dec 09 '12 at 11:19