2
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ 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". $\endgroup$ Commented Feb 21, 2013 at 19:10
  • $\begingroup$ What is the bit rate / bit duration? $\endgroup$ Commented Feb 21, 2013 at 22:09
  • $\begingroup$ 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. $\endgroup$ Commented Feb 21, 2013 at 22:53
  • $\begingroup$ Try using fftshift(). I think it would help your issue. $\endgroup$ Commented Feb 2, 2014 at 7:13

1 Answer 1

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Feb 22, 2013 at 7:55
  • 1
    $\begingroup$ A common way to avoid linspace with Octave/Matlab is using t = Ts * (0:N-1) $\endgroup$ Commented Feb 22, 2013 at 15:56
  • 1
    $\begingroup$ NumPy's linspace is easier: you just say endpoint=False :) $\endgroup$ Commented Feb 26, 2013 at 22:21

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.