Hi,
I’ve been bashing my head for a couple of days on what should be a
trivial
process.
I would really appreciate it if somebody could please come to my rescue.
I’d like to digitise a signal, intented to be off air FM signal, but for
the
moment, i’m capturing a 50KHz sinusiod and doing the reconstruction in
MATLAB… to test that im capturing the samples correctly.
on plotting the FFT in MATLAB, it’s clear that something’s a miss…I
don’t
know where 
i’ve been at this for ages!!!
here’s my code: I really hope somebody can help, I know this should be a
simple matter. The Matlab code follows below the Python stuff.
=============================
#!/usr/bin/env python
“”"
Read samples from the USRP and write to file formatted as binary
outputs single precision complex float values or complex short values
(interleaved 16 bit signed short integers).
“”"
from gnuradio import gr, eng_notation
from gnuradio import audio
from gnuradio import usrp
from gnuradio.eng_option import eng_option
from optparse import OptionParser
class my_graph(gr.flow_graph):
def __init__(self):
gr.flow_graph.__init__(self)
usage="%prog: [options] output_filename"
parser = OptionParser(option_class=eng_option, usage=usage)
parser.add_option("-R", "--rx-subdev-spec", type="subdev",
default=(0, 0),
help=“select USRP Rx side A or B (default=A)”)
parser.add_option("-d", “–decim”, type=“int”, default=8,
help=“set fgpa decimation rate to DECIM
[default=%default]”)
parser.add_option("-f", “–freq”, type=“eng_float”,
default=50e3,
help=“set frequency to FREQ”, metavar=“FREQ”)
parser.add_option("-g", “–gain”, type=“eng_float”,
default=None,
help=“set gain in dB (default is midpoint)”)
parser.add_option("-N", “–nsamples”, type=“eng_float”,
default=2000,
help=“number of samples to collect
[default=+inf]”)
(options, args) = parser.parse_args ()
if len(args) != 1:
parser.print_help()
raise SystemExit, 1
filename = args[0]
if options.freq is None:
parser.print_help()
sys.stderr.write('You must specify the frequency with -f
FREQ\n’);
raise SystemExit, 1
# build the graph
self.u = usrp.source_c(decim_rate=options.decim)
self.dst = gr.file_sink(gr.sizeof_gr_complex, filename)
self.head = gr.head(gr.sizeof_gr_complex, int(options.nsamples))
self.connect(self.u, self.head, self.dst)
rx_subdev_spec = usrp.pick_rx_subdevice(self.u)
self.u.set_mux(usrp.determine_rx_mux_value(self.u,
options.rx_subdev_spec))
# determine the daughterboard subdevice we're using
self.subdev = usrp.selected_subdev(self.u,
options.rx_subdev_spec)
print “Using RX d’board %s” % (self.subdev.side_and_name(),)
input_rate = self.u.adc_freq() / self.u.decim_rate()
print “USB sample rate %s” %
(eng_notation.num_to_str(input_rate))
print “Freq is set to: %s” % (options.freq)
if options.gain is None:
# if no gain was specified, use the mid-point in dB
g = self.subdev.gain_range()
options.gain = float(g[0]+g[1])/2
self.subdev.set_gain(options.gain)
r = self.u.set_rx_freq (0, options.freq) #self.u.tune(0,
self.subdev,
options.freq)
if not r:
sys.stderr.write(‘Failed to set frequency\n’)
raise SystemExit, 1
if name == ‘main’:
try:
my_graph().run()
except KeyboardInterrupt:
pass
MATLAB
ms = 1e-3;
kHz = 1e3;
count = 2000;
decim_rate = 8;
Fsamp = 64e6;
Fsamp_real = Fsamp./decim_rate;
%set the time axis
dt = 1./Fsamp_real;
T_end = count./Fsamp_real;
t = 0:dt:T_end-dt;
%t = dt:dt:T_end;
df = fopen(‘signal_samples.dat’);
y = fread (df,[2, count], ‘float’);
fclose(df);
%plot I and Q
figure (1)
subplot(211)
plot(t/ms,y(1,:));
grid on;
xlabel(‘ms’)
title (‘I data’)
subplot(212)
plot(t/ms,y(2,:));
grid on;
xlabel(‘ms’)
title (‘Q data’)
%remove the artefact
z = y(:,601:end);
t0 = t(1:end-600); %new time!
%plot I and Q again
figure (2)
subplot(211)
plot(t0/ms,z(1,:));
grid on;
xlabel(‘ms’)
title (‘I data’)
subplot(212)
plot(t0/ms,z(2,:));
grid on;
xlabel(‘ms’)
title (‘Q data’)
%compute fft
z_spectrum = fftshift(fft(z));
%compute frequency axis
df = 1/(t0(end)-t0(1));
freq_axis = -Fsamp_real/2 :df :Fsamp_real/2;
%FFT module
figure (3)
subplot(211)
plot(freq_axis/kHz,abs(z_spectrum(1,:)))
xlabel(‘kHz’)
title (‘I data spectrum’)
subplot(212)
plot(freq_axis/kHz,abs(z_spectrum(2,:)))
xlabel(‘kHz’)
title (‘Q data spectrum’)
%FFT phase
figure (4);hold on;
plot(freq_axis/kHz,(angle(z_spectrum(1,:))*180/pi))
plot(freq_axis/kHz,(angle(z_spectrum(2,:))*180/pi),‘r’)
xlabel(‘kHz’)
ylabel(‘deg’)
title (‘I and Q phases’)
==========================================
i’ve been at this for ages!!!
i’ve been at this for ages!!!