2

I'm writing a function that takes in a time-varying signal and returns the FFT. It follows the Matlab documentation:

It works well for really simple sums of sinusoids.

For other simulated signals, like step functions, I end up with something fubar, that oscillates in a sawtooth-like fashion. Here is a minimal example using a heaviside step function, that plots the imaginary component of the FFT weight:

%Create and plot sign function 
Fs = 1000;
dt = 1/Fs;
tvals = -1: dt: 1-dt;
yvals = sign(tvals); %don't use heaviside() it is very slow

subplot(2,1,1);
    plot(tvals, yvals, 'k', 'LineWidth', 2); hold on; 
    axis([-1 1 -1.1 1.1]);

%Calculate center-shifted FFT and plot imaginary part
fftInitial = fft(yvals);
n = length(fftInitial);
dF = Fs/n;
frequenciesShifted = -Fs/2: dF: Fs/2-dF;  %zero-centered frequency range
fftShifted = fftshift(fftInitial/length(yvals));

subplot(2,1,2);
plot(frequenciesShifted, imag(fftShifted), 'b', 'LineWidth', 2);hold on
xlim([-8 8])

Here is the resulting plot: enter image description here

Note the imaginary solution is known to be 2/jw, or j(-2/w).

Note in addition to the sawtooth, which is my primary concern, the envelope of these weights doesn't even seem to follow that. It actually seems flipped around the origin. Not sure what combination of things I've done wrong here.

Based on some helpful feedback, people pointed out this question:
Analytical Fourier transform vs FFT of functions in Matlab

In particular, not including a zero in the time array could cause problems, which was a main concern there. The time array in my code includes a zero, so it doesn't seem like a duplicate. I have pushed my zero to the front of the array by time-shifting the step (though frankly I have done lots of fft() without doing this and they have all looked fine, so I don't imagine that is the problem, but I just want to remove that as an issue for now). So we end up with:

Fs = 1000;
dt = 1/Fs;
tvals = 0: dt: 2-dt;
yvals = sign(tvals-1); %don't use heaviside() it is very slow
zeroInd = find(yvals == 1, 1, 'first');
yvals(zeroInd-1) =0;
%Calculate center-shifted FFT and plot imaginary part
fftInitial = fft(yvals);
n = length(fftInitial);
dF = Fs/n;
frequenciesShifted = -Fs/2: dF: Fs/2-dF;  %zero-centered frequency range
fftShifted = fftshift(fftInitial/length(yvals));

%Plot stuff
subplot(2,1,1);
plot(tvals, yvals, 'k', 'LineWidth', 2); hold on; 
axis([0 2 -1.1 1.1]);  
subplot(2,1,2);
plot(frequenciesShifted, imag(fftShifted), 'b', 'LineWidth', 2);hold on
xlim([-8 8])
grid on;

I still get the same zig-zagging function with the same incorrect envelope. So, while I do see that my question and that question are closely related, I am not sure they are duplicates. And I would really like to be able to draw a halfway-reasonable looking plot of the values of these components (I do work with some physiological signals that are basically step functions (they move very fast compared to my measuring instruments, so this is not just an academic exercise for me).

eric
  • 7,142
  • 12
  • 72
  • 138
  • 1
    The duplicate question uses a different function, but has exactly the same issue you are facing: the FFT assumes that the origin (t=0) is the first element in the input array. You are computing the FFT of a shifted step function. – Cris Luengo Sep 13 '18 at 17:46
  • That question is not really a duplicate: for one, I have a zero in my time representation, and that was the problem in their post. Second, my main problem is the oscillatory behavior of the FT, not the mirror reflection, as indicated in my title. But I will edit my Q to mention that post. – eric Sep 13 '18 at 18:06
  • @CrisLuengo Cris I read it:...As you wrote there "The FFT expects the origin to be at the first (leftmost) sample. This is what ifftshift is for:" So what I did in my revised question was just make 0 my first time point. I thought that would obviate the need for iffshift. Will look more closely at this now. – eric Sep 13 '18 at 18:48
  • @CrisLuengo why does this extremely simple solution work for a unit pulse function (without needing to muck about with time and such), with a gorgeous plot, but when I do the exact same thing with a single step function, it is a disaster? https://stackoverflow.com/questions/19717779/i-did-a-step-function-fft-on-matlab-but-only-getting-one-side-of-the-spectrum-0 That is, why is the answer there so simple and easy. I'm just not seeing the difference (especially given that a pulse is just a sum of two signums!). – eric Sep 13 '18 at 19:25
  • I was wrong, you have that issue but also a second one. Sorry. See my answer. The question with the unit pulse you link here plots the magnitude of the Fourier transform, so shifts in the signal cannot be seen. – Cris Luengo Sep 13 '18 at 19:29

1 Answers1

3

There are two problems with your code:

  1. The DFT (the FFT algorithm computes the DFT) takes the origin to be at the leftmost bin. Creating yvals such that the origin is in the middle causes the output Frequency spectrum to be that of a signal shifted by half its length. This leads to very high frequency oscillations. The fix is to use ifftshift on the input data before calling fft. See more in this other question

  2. The DFT assumes (can be interpreted as assuming) that the input signal is periodic. This leads to a second large jump. Basically, your signal looks like a shifted box function, and hence your transform looks like a sinc function with modified phase. The solution is to apply a windowing function to your input before calling fft. See for example this other question.

Modify your code as follows:

yvals = sign(tvals);
yvals = yvals .* hanning(numel(yvals), 'periodic').'; % Apply windowing function

% ...

fftInitial = fft(ifftshift(yvals)); % Shift signal before calling FFT

This is the output your code now gives:

graph

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • 1
    So basically I picked one of the worst possible functions to start building my intuitions about FFT. :) – eric Sep 13 '18 at 19:59