# 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

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
NumPy's linspace is easier: you just say endpoint=False :) – endolith Feb 26 '13 at 22:21