spec_analyze.m

% Spectrum Analysis with RTL-SDR Radio

FILEMODE = 1; % If this flag is set to 0, uses "real" data from RTL-SDR. If set to 1, uses captured data in lab2.dat instead.

% To complete the lab, you shouldn't change these parameters (unless you want to mess
% around and explore beyond the lab instructions... that's encouraged!
% If you do so, include any observations in your lab report!)
fc = 910.0e6; % Center frequency (Hz)
Ts = 1/2e6;   % Samples per second
FrameLength = 256*20;  % number of samples to "grab" each time through loop
simLength = 5000;      % total number of frames to grab (determine total sim time)
                    % Note: total simulation time = simLength * FrameLength * Ts

% Create receiver object
if FILEMODE == 0
    hSDRrRx = comm.SDRRTLReceiver(...
        'CenterFrequency', fc, ...
        'EnableTunerAGC',  true, ...
        'SampleRate',      round(1/Ts), ...
        'SamplesPerFrame', FrameLength, ...
        'OutputDataType',  'double');
else % unless we're in file mode, then open file and set simLength appropriately
    simLength = 1000;
    fid=fopen('lab2.dat');
end

% create spectrum analyzer object
hSpectrum = dsp.SpectrumAnalyzer(...
    'Name',             'Baseband Spectrum',...
    'Title',            'Baseband Spectrum', ...
    'SpectrumType',     'Power density',...
    'FrequencySpan',    'Full', ...
    'SampleRate',       round(1/Ts), ...
    'YLimits',          [-80,10],...
    'SpectralAverages', 50, ...
    'FrequencySpan',    'Start and stop frequencies', ...
    'StartFrequency',   -round(1/Ts/2), ...
    'StopFrequency',    round(1/Ts/2),...
    'Position',         figposition([50 30 30 40]));

theta1 = zeros(1,5120);
theta2 = zeros(1,5120);
mu1=5; mu2=.005;
f0=800e3/2;
 
% Main loop to grab samples
for count = 1 : simLength
    if FILEMODE == 0
        [data, ~] = step(hSDRrRx);       % grab complex (i.e. quadrature) samples from RTL-SDR
        data = real(data - mean(data));  % remove DC component, and only keep real portion
    else % grab data from file instead
        data=fread(fid, 5120, 'double');
        
        %pause(0.003); % slow things down a little to mimic real-time
    end
    squaredData=data.^2;
    
    fcenter=800e3;
    fspan=.00008*910e6*2;
    fl=100; ff=[0  fcenter-fspan-25e3 fcenter-72.8e3 fcenter+72.8e3 fcenter+fspan+25e3 1e6]/1e6; 
    fa=[0 0 1 1 0 0];
    b = firpm(fl, ff, fa);
    if exist('lastData') 
        tempData = filter(b, 1, [lastData; squaredData]);
        filteredData = tempData(end-5120:end);
    else 
        filteredData = filter(b, 1, data);
    end
        
    lastData = data;
    lent=5120; carest=zeros(1,lent);
    t=(count-1)*5120*Ts:Ts:5120*count*Ts;
    for k=1:lent-1                      % combine top PLL theta1
        theta1(k+1)=theta1(k)-mu1*filteredData(k)*sin(4*pi*600e3*t(k)+2*theta1(k));           
        theta2(k+1)=theta2(k)-mu2*filteredData(k)*sin(4*pi*600e3*t(k)+2*theta1(k)+2*theta2(k));  
        % with bottom PLL theta2 to form estimate of preprocessed signal
        carest(k)=cos(2*pi*600e3*t(k)+theta1(k)+theta2(k));
        theta1(1) = theta1(end);
        theta2(1) = theta2(end);
    end
    fl=100; ff=[0 .2 .3 1]; fa=[1 1 0 0];
    b = firpm(fl, ff, fa);
    %output = filter(b, 1, data.*carest');
    step(hSpectrum, data.*carest');           % update spectrum analyzer display
end
    
if FILEMODE == 0  % close RTL-SDR object
    release(hSDRrRx);
else % close file
    fclose(fid);
end

% Release all system objects
release(hSpectrum);

subplot(2, 1, 1);
plot(theta1(theta1 ~= 0))
subplot(2, 1, 2);
%hold on
plot(theta2(theta2 ~= 0))
hold off