FFT in MATLAB is shifted in frequency

Question

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.

Answer 1

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.