Converting Scales after FFT

I was hoping someone could help me understand exactly how Matlab's FFT function deals with scaling axes. I am developing some Diffraction code that will simulate Fourier optics but I am hoping to understand how to interpret my results better with a simple example first.

Analytically the Fourier Transform of a f(x) = cos(bx) wave goes to F(k) = dirac(k-b)+dirac(k+b). My code to test this in Matlab is:

steps = 2^8; lim = 4*pi;

x=linspace(-lim, lim,steps);

b = 1;
func = cos(b*x);
Func = fftshift(fft(func));

which when plotted looks like this:

What I am trying to understand is how to rescale my k (spatial frequency) axis so the locations of the delta functions match the analytic result, which in this case should b at k = +/- 1.

This is my first time posting to these forums. Please inform and forgive me if I break any conventions. Any advice you could offer would be greatly appreciated!

Answer 1

[Star Strider]

The current R2015b documentation for fft is confusing (at least in my opinion). The R2015a documentation is clearer: fft, but still contains a scaling error w.r.t. plotting a one-sided fft, leaving out the 2 factor.

Calculating the fft as:

Func = fft(func)*2/length(func);

will produce the correct magnitudes.

[James Kirk]

Thank you very much for your answer. I see how the code you suggested normalises the magnitude of the peaks. I am still confused however as to how to interpret the separation of the peaks. Are you able to explain how I could construct a new x/k scale that would line up with the analytic result?

This also applies to my other question. Where I am confused as to why he number of steps alters the width of the transformed function.

Thank you for your help so far!

[Star Strider]

There is actually only one peak. The frequency axis is imaginary and so the plot of the full fft plots the complex conjugate frequencies, ±jω. For clarity, and since the -jω plot is the mirror-image of the +jω plot, usually the +jω section is the only one plotted.

The full code to do that would be:

steps = 2^8; lim = 4*pi;
x=linspace(-lim, lim,steps);

b = 1;
func = cos(b*x);

Ts = mean(diff(x));                                 % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency
Fn = Fs/2;                                          % Nyquist Frequency
Func = fft(func)*2/steps;                           % Normalised FFT
Fv = linspace(0, 1, fix(steps/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                                  % Index Vector (For Plotting)

figure(1)
plot(Fv, abs(Func(Iv)))
grid
axis([0  1    ylim])                                % Limit Axis Displayed
xlabel('Frequency (Hz)')
ylabel('Amplitude')