# FFT in MATLAB is shifted in frequency

I have a signal in MATLAB on which I'd like to perform an FFT. My signal is stored in s, and I use the code below (inspired by the MATLAB help):

L = length(s);
nfft = 2^nextpow2(L);
S = fft(s,nfft)/L;
fftf = 1/Ts/2*linspace(0,1,nfft/2+1);
ffts = (2*abs(S(1:nfft/2+1)));


s is a GMSK-modulated bit sequence, that is, it varies between $f_c-2400$ Hz and $f_c+2400$ Hz in my case when transmitting 0 or 1, respectively. $f_c$ is set to 100kHz. For a long input, say, 250 bits worth of 1's and 250 bits worth of 0's, I get what I expect, see the first image below.

If I instead choose a low number of bits, say 10 1 bits followed by 10 0 bits, I get as expected, but it is shifted down to ~90kHz instead of centered at 100kHz. This is something I can't quite understand - it seems changing the sample rate and length of the FFT changes absolutely nothing.

Can anyone explain to me why? Thanks in advance!

Long data: Short data: The code used to generate the signal and FFT:

%% Configuration
clear; clf;

DataRate = 9600;          % 9600 kbps for AIS
N = 100;                  % Samples per bit
Tb = 1/DataRate;          % Bit period
Ts = Tb/N;                % Sampling period
BT = 0.5;                 % AIS spec, time-bandwidth product
Ftrans = 100e3;           % Frequency of "transmitted" signal

num = 200;
Bits = zeros(1,num)+1;
Bits = [Bits zeros(1,num)-1];
clear num;

%% Modulation

% Prep a time axis from -2Tb to 2Tb
t_g = -2*Tb:Ts:2*Tb;

% Gaussian response to rectangular pulse [Haykin4th, p. 397]
x = pi*sqrt(2/log(2))*BT;
gr = 1/2*(erfc(x*(t_g/Tb-1/2))-erfc(x*(t_g/Tb+1/2)));

% Truncate to 3Tb, pulse centered at 1.5Tb
gr = gr(0.5*N+2:3.5*N+1);

% Normalize
% when integrating, we want to end at 0.5 (phase changes by 0.5pi)
% so, we want sum(y)=0.5 -> normalize by sum(y) and divide by two.
gr = (gr)./(2*sum(gr));

% Generate the Gaussian filtered pulse train by centering a "Gaussian
% rectangle" on each bit, and adding inter-symbol interference
f = zeros(1,(length(Bits)+2)*N);
for n = 1:length(Bits)
f((n-1)*N+1:(n+2)*N) = f((n-1)*N+1:(n+2)*N) + Bits(n).*gr;
end

% Since gr corresponds to changing the phase 0.5, multiplying by pi and
% integrating gives the resulting phase.
theta = pi*cumsum(f);

% Prep I,Q
I = cos(theta);
Q = sin(theta);

% Transmitted signal, shifted to ftrans
t = linspace(0,length(Bits)*N,length(I))*Ts;
s = -sin(2*pi*Ftrans.*t).*Q+cos(2*pi*Ftrans.*t).*I;

%% FFT

L = length(s);
% faster w/ a pow2 length, signal padded with zeros
nfft = 2^nextpow2(L);
% do the nfft-point fft and normalize
S = fft(s,nfft)/L;
% x-axis from 0 to fs/2, nfft/2+1 points
fftf = 1/Ts/2*linspace(0,1,nfft/2+1);
% only plotting the first half since its mirrored, thus 1:nfft/2+1
% why multiplied with 2?
ffts = (2*abs(S(1:nfft/2+1)));

%% Plotting

% FFT PLOT
plot(fftf/1e3,ffts);
title('FFT of transmitted signal S(f)');
set(gca,'xlim',[Ftrans/1e3-20 Ftrans/1e3+20]);
ylabel('|S(f)|');
xlabel('Frequency [kHz]');
grid;


Adjusting the sample frequency by changing N seems to have no effect - but changing num from e.g. 10 to 100 (changing the number of bits) clearly shifts the plotted spectrum closer to 100kHz.

• Can you post your signal generation code? The two tones don't seem to be well centered about 100 kHz, which makes me think there might be a subtle bug in the signal generation code. Also, the tones around 90 kHz are off in the "other direction". – Dave C Feb 21 '13 at 19:10
• What is the bit rate / bit duration? – user2718 Feb 21 '13 at 22:09
• Thanks for the replies. I added the code above - I've inspected the f, theta, I, Q and s signals quite a bit and they seem pretty good to me. Please, let me know if you think otherwise :) The bit rate is 9600bps. – Tausen Feb 21 '13 at 22:53
• Try using fftshift(). I think it would help your issue. – sundar Feb 2 '14 at 7:13

## 1 Answer

Unfortunately I don't have Matlab now but I'm thinking your problem has to do with linspace. Remember linspace gives the actual first and last values always. When you're doing the FFT binning I don't believe you want the boundaries to be pinned that way. I'd try using the colon operator instead and appropriate starting and delta values.

EDIT: If you look at the definition of the DFT it is: $$X(k)=\sum_{n=0}^{N−1}x(n)e^{j\frac{2πkn}{N}}$$

This shows that the first DFT output at $k=0$ is at baseband. Every subsequent frequency is at an interval of $F_s/N$. So if you had 10 points the bin centers would be from $0...(\frac{N-1}{N}F_s=0.9F_s)$. For 1000 points it would span $0...0.999F_s$. Linspace, on the other hand always leaves the first and last points exact. You need the start and spacing exact, or to adjust the ending point with linspace.

• Fantastic! You are absolutely right - using this time axis solved the problem: t = 0:Ts:(length(Bits)+2)*Tb-Ts; I'm not sure I understand why the linspace axis doesn't work, though - guess its back to the books for me. – Tausen Feb 22 '13 at 7:55
• A common way to avoid linspace with Octave/Matlab is using t = Ts * (0:N-1) – Dave C Feb 22 '13 at 15:56
• NumPy's linspace is easier: you just say endpoint=False :) – endolith Feb 26 '13 at 22:21