Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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:

http://f.cl.ly/items/0R1i0x0m46221A230m25/fft1.png

Short data:

http://f.cl.ly/items/2a320I2l2r3M1i2o2N11/fft1.png

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.

share|improve this question
    
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 at 7:13
add comment

1 Answer

up vote 6 down vote accepted

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.

share|improve this answer
1  
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
1  
A common way to avoid linspace with Octave/Matlab is using t = Ts * (0:N-1) –  Dave C Feb 22 '13 at 15:56
1  
NumPy's linspace is easier: you just say endpoint=False :) –  endolith Feb 26 '13 at 22:21
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.