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.