Adding a Controlled Amount of Noise to a Noise-Free Signal

Rick Lyons January 3, 201221 comments Coded in Matlab
function [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
%   SNR_Set(x, Desired_SNR_dB) returns the real-valued 
%   input 'Signal' contaminated with normally-distributed, 
%   zero-mean, random noise.  The signal-to-noise ratio   
%   (SNR in dB) of the output 'Noisy_Signal' signal is  
%   controlled by the input 'Desired_SNR_dB' argument measured 
%   in dB.

%   Example:
% 	Npts = 128;                  % Number of time samples
% 	n = 0:Npts-1;                % Time-domain index
% 	Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% 	Desired_SNR_dB = 3;      % Set SNR of output 'Noisy_Signal' to +3 dB
% 	[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
%   Author: Richard Lyons [December 2011]
%******************************************

Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one

Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
	%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));

K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10);  % Scale factor

New_Noise = sqrt(K)*Noise; % Change Noise level
	%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
	%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))

Noisy_Signal = Signal + New_Noise;

'SNR_Set()' Function Test Code:
%  Filename SNR_Set_test.m
%
%  Tests the 'SNR_Set()" function.  Adds a predefined 
%  amount of random noise to a noise-free signal such that  
%  the noisy signal has a desired signal-to-noise ratio (SNR).
%
%  Author: Richard Lyons [December 2011]

clear, clc

% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5;  % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal

Desired_SNR_dB = 3 % Set desired SNR in dB

[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal

% Plot original and 'noisy' signals
    figure(1)
    subplot(2,1,1)
    plot(n, Signal, '-bo', 'markersize', 2)
    title('Original Signal')
    grid on, zoom on
    subplot(2,1,2)
    plot(n, Noisy_Signal, '-bo', 'markersize', 2)
    title('Noisy Signal')
    xlabel('Time-samples')
    grid on, zoom on

% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
    figure(2)
    plot(Spec_Mag, '-bo', 'markersize', 2)
    title('Spec Mag of Noisy Signal')
    xlabel('Freq-samples'), ylabel('Linear')
    grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)

Narrow-band moving-average decimator, one addition/sample

Markus Nentwig December 31, 2011 Coded in Matlab
% *************************************************
% Moving average decimator 
%
% A decimator for narrow-band signals (~5 % or less bandwidth occupation) 
% can be implemented as follows:
%
% #define DECIM (100)
% double acc = 0.0;
% while (1){
%    int ix;
%    for(ix = DECIM; ix > 0; --ix){
%       acc += getInputSample();
%    } /* for */
%    writeOutputSample(acc / (double)DECIM);
%    acc = 0.0;
% } /* while */
%
% It is conceptually identical to a moving average filter
% http://www.dspguide.com/ch15/4.htm combined with a decimator
%
% Note that the "moving" average jumps ahead in steps of the decimation
% factor. Intermediate output is decimated away, allowing for a very efficient 
% implementation.
% This program calculates the frequency response and alias response,
% based on the decimation factor and bandwidth of the processed signal.
% *************************************************    
function eval_design()
    decimationFactor = 100;
    rateIn_Hz = 48000;
    noDecim = false;
    
    % create illustration with sinc response
    %decimationFactor = 4; noDecim = true; 
    
    % *************************************************
    % signal source: Bandlimited test pulse
    % Does not contain energy in frequency ranges that 
    % cause aliasing
    % *************************************************
    s = zeros(1, 10000 * decimationFactor);    
    fb = FFT_frequencyBasis(numel(s), rateIn_Hz);
    
    % assign energy to frequency bins
    if noDecim
        sPass = ones(size(s));
    else
        sPass = s; sPass(find(abs(fb) < rateIn_Hz / decimationFactor / 2)) = 1;
    end
    sAlias = s; sAlias(find(abs(fb) >= rateIn_Hz / decimationFactor / 2)) = 1;
    
    % convert to time domain
    sPass = fftshift(real(ifft(sPass)));
    sAlias = fftshift(real(ifft(sAlias)));

    % *************************************************    
    % plot spectrum at input
    % *************************************************    
    pPass = {};
    pPass = addPlot(pPass, sPass, rateIn_Hz, 'k', 5, ...
                    'input (passband response)');
    pAlias = {};
    pAlias = addPlot(pAlias, sAlias, rateIn_Hz, 'k', 5, ...
                     'input (alias response)');    
        
    % *************************************************    
    % impulse response
    % *************************************************    
    h = zeros(size(s));
    h (1:decimationFactor) = 1; 
    if noDecim
        h = h / decimationFactor;
        decimationFactor = 1;
    end
    
    % cyclic convolution between signal and impulse response
    sPass = real(ifft(fft(sPass) .* fft(h)));
    sAlias = real(ifft(fft(sAlias) .* fft(h)));
    
    % decimation
    sPass = sPass(decimationFactor:decimationFactor:end);
    sAlias = sAlias(decimationFactor:decimationFactor:end);
    rateOut_Hz = rateIn_Hz / decimationFactor;
    
    % *************************************************    
    % plot spectrum
    % *************************************************    
    pPass = addPlot(pPass, sPass, rateOut_Hz, 'b', 3, ...
                    'decimated (passband response)');
    figure(1); clf; grid on; hold on;
    doplot(pPass, sprintf('passband frequency response over input rate, decim=%i', decimationFactor));
    
    pAlias = addPlot(pAlias, sAlias, rateOut_Hz, 'b', 3, ...
                     'decimated (alias response)');
    figure(2); clf; grid on; hold on;
    doplot(pAlias, sprintf('alias frequency response over input rate, decim=%i', decimationFactor));

    % *************************************************    
    % plot passband ripple
    % *************************************************    
    fb = FFT_frequencyBasis(numel(sPass), 1);
    fr = 20*log10(abs(fft(sPass) + eps));
    ix = find(fb > 0);

    figure(3); clf;
    h = semilogx(fb(ix), fr(ix), 'k');
    set(h, 'lineWidth', 3);
    ylim([-3, 0]);
    title(sprintf('passband gain over output rate, decim=%i', decimationFactor));
    xlabel('frequency relative to output rate');
    ylabel('dB'); grid on;
    
    % *************************************************    
    % plot alias response
    % *************************************************    
    fb = FFT_frequencyBasis(numel(sAlias), 1);
    fr = 20*log10(abs(fft(sAlias) + eps));
    ix = find(fb > 0);
    
    figure(4); clf;
    h = semilogx(fb(ix), fr(ix), 'k');
    set(h, 'lineWidth', 3);
    %    ylim([-80, -20]);
    title(sprintf('alias response over output rate, decim=%i', decimationFactor));
    xlabel('frequency relative to output rate');
    ylabel('dB'); grid on;    
end

% ************************************
% put frequency response plot data into p
% ************************************
function p = addPlot(p, s, rate, plotstyle, linewidth, legtext)
    p{end+1} = struct('sig', s, 'rate', rate, 'plotstyle', plotstyle, 'linewidth', linewidth, 'legtext', legtext);
end

% ************************************
% helper function, plot data in p
% ************************************
function doplot(p, t)
    leg = {};
    for ix = 1:numel(p)
        pp = p{ix};        
        fb = FFT_frequencyBasis(numel(pp.sig), pp.rate);
        fr = 20*log10(abs(fft(pp.sig) + eps));
        h = plot(fftshift(fb), fftshift(fr), pp.plotstyle);
        set(h, 'lineWidth', pp.linewidth);
        xlabel('f / Hz');
        ylabel('dB');
        leg{end+1} = pp.legtext;
    end
    legend(leg);
    title(t);
end

% ************************************
% calculates the frequency that corresponds to
% each FFT bin (negative, zero, positive)
% ************************************
function fb_Hz = FFT_frequencyBasis(n, rate_Hz)
    fb = 0:(n - 1);
    fb = fb + floor(n / 2);
    fb = mod(fb, n);
    fb = fb - floor(n / 2);
    fb = fb / n; % now [0..0.5[, [-0.5..0[
    fb_Hz = fb * rate_Hz;
end

Alias/error simulation of interpolating RRC filter

Markus Nentwig December 25, 2011 Coded in Matlab
% *************************************************
% alias- and in-channel error analysis for root-raised
% cosine filter with upsampler FIR cascade
% Markus Nentwig, 25.12.2011
% 
% * plots the aliases at the output of each stage
% * plots the error spectrum (deviation from ideal RRC-
%   response)
% *************************************************    
function eval_RRC_resampler()
    1;
    % variant 1
    % conventional RRC filter and FIR resampler
    smode = 'evalConventional'; 
    
    % export resampler frequency response to design equalizing RRC filter
    %smode = 'evalIdeal';
    
    % variant 2
    % equalizing RRC filter and FIR resampler
    % smode = 'evalEqualized';
    
    % *************************************************
    % load impulse responses
    % *************************************************    
    switch smode
      case 'evalConventional'
        % conventionally designed RRC filter
        h0 = load('RRC.dat');
      
      case 'evalIdeal'
        h0 = 1;
      
      case 'evalEqualized'
        % alternative RRC design that equalizes the known frequency 
        % response of the resampler
        h0 = load('RRC_equalized.dat');
      
      otherwise assert(false);
    end
        
    h1 = load('ip1.dat');
    h2 = load('ip2.dat');
    h3 = load('ip3.dat');
    h4 = load('ip4.dat');

    % *************************************************
    % --- signal source ---
    % *************************************************
    n = 10000; % test signal, number of symbol lengths 
    rate = 1;    
    s = zeros(1, n);
    s(1) = 1;

    p = {};
    p = addPlot(p, s, rate, 'k', 5, 'sym stream r=1');

    % *************************************************
    % --- upsampling RRC filter ---
    % *************************************************
    rate = rate * 2;
    s = upsample(s, 2); % insert one zero after every sample
    p = addPlot(p, s, rate, 'k', 2, 'sym stream r=2');
    s = filter(h0, [1], s);
    p = addPlot(p, s, rate, 'b', 3, 'sym stream r=2, filtered');
    p = addErrPlot(p, s, rate, 'error');

    figure(1); clf; grid on; hold on;
    doplot(p, 'interpolating pulse shaping filter');
    ylim([-70, 2]);
    p = {};

    % *************************************************
    % --- first interpolator by 2 ---
    % *************************************************
    rate = rate * 2;
    s = upsample(s, 2); % insert one zero after every sample
    p = addPlot(p, s, rate, 'k', 3, 'interpolator 1 input');
    s = filter(h1, [1], s);
    p = addPlot(p, s, rate, 'b', 3, 'interpolator 1 output');
    p = addErrPlot(p, s, rate, 'error');

    figure(2); clf; grid on; hold on;
    doplot(p, 'first interpolator by 2');
    ylim([-70, 2]);
    p = {};
    
    % *************************************************
    % --- second interpolator by 2 ---
    % *************************************************
    rate = rate * 2;
    s = upsample(s, 2); % insert one zero after every sample
    p = addPlot(p, s, rate, 'k', 3, 'interpolator 2 input');
    s = filter(h2, [1], s);
    p = addPlot(p, s, rate, 'b', 3, 'interpolator 2 output');
    p = addErrPlot(p, s, rate, 'error');

    figure(3); clf; grid on; hold on;
    doplot(p, 'second interpolator by 2');
    ylim([-70, 2]);
    p = {};

    % *************************************************
    % --- third interpolator by 2 ---
    % *************************************************
    rate = rate * 2;
    s = upsample(s, 2); % insert one zero after every sample
    p = addPlot(p, s, rate, 'k', 3, 'interpolator 3 input');
    s = filter(h3, [1], s);
    p = addPlot(p, s, rate, 'b', 3, 'interpolator 3 output');
    p = addErrPlot(p, s, rate, 'error');

    figure(4); clf; grid on; hold on;
    doplot(p, 'third interpolator by 2');
    ylim([-70, 2]);
    p = {};
    
    % *************************************************
    % --- fourth interpolator by 4 ---
    % *************************************************
    rate = rate * 4;
    s = upsample(s, 4); % insert three zeros after every sample
    p = addPlot(p, s, rate, 'k', 3, 'interpolator 4 input');
    s = filter(h4, [1], s);
    p = addPlot(p, s, rate, 'b', 3, 'final output');
    p = addErrPlot(p, s, rate, 'error at output');
    
    figure(5); clf; grid on; hold on;
    doplot(p, 'fourth interpolator by 4');
    ylim([-70, 2]);    
    
    figure(334);
    stem(real(s(1:10000)));
    
    % *************************************************
    % export resampler frequency response
    % *************************************************
    switch smode
      case 'evalIdeal'
        exportFrequencyResponse(s, rate, 'interpolatorFrequencyResponse.mat');
    end    
end

% ************************************
% put frequency response plot data into p
% ************************************
function p = addPlot(p, s, rate, plotstyle, linewidth, legtext)
    p{end+1} = struct('sig', s, 'rate', rate, 'plotstyle', plotstyle, 'linewidth', linewidth, 'legtext', legtext);
end

% ************************************
% determine the error spectrum, compared to an ideal filter (RRC)
% and add a plot to p
% ************************************
function p = addErrPlot(p, s, rate, legtext)

    ref = RRC_impulseResponse(numel(s), rate);
    % refB is scaled and shifted (sub-sample resolution) replica of ref
    % that minimizes the energy in (s - refB)
    [coeff, refB, deltaN] = fitSignal_FFT(s, ref);
    err = s - refB;
    err = brickwallFilter(err, rate, 1.15); % 1+alpha

    % signal is divided into three parts:
    % - A) wanted in-channel energy (correlated with ref)
    % - B) unwanted in-channel energy (uncorrelated with ref)
    % - C) unwanted out-of-channel energy (aliases)
    % the error vector magnitude is B) relative to A) 
    energySig = refB * refB';
    energyErr = err * err';
    EVM_dB = 10*log10(energyErr / energySig);
    legtext = sprintf('%s; EVM=%1.2f dB', legtext, EVM_dB);
    
    p{end+1} = struct('sig', err, 'rate', rate, 'plotstyle', 'r', 'linewidth', 3, 'legtext', legtext);
end

% ************************************
% helper function, plot data in p
% ************************************
function doplot(p, t)
    leg = {};
    for ix = 1:numel(p)
        pp = p{ix};        
        fb = FFT_frequencyBasis(numel(pp.sig), pp.rate);
        fr = 20*log10(abs(fft(pp.sig) + eps));
        h = plot(fftshift(fb), fftshift(fr), pp.plotstyle);
        set(h, 'lineWidth', pp.linewidth);
        xlabel('frequency, normalized to symbol rate');
        ylabel('dB');
        leg{end+1} = pp.legtext;
    end
    legend(leg);
    title(t);
end

% ************************************
% ideal RRC filter (impulse response is as 
% long as test signal)
% ************************************
function ir = RRC_impulseResponse(n, rate)
    alpha = 0.15;
    fb = FFT_frequencyBasis(n, rate);
    % bandwidth is 1
    c = abs(fb / 0.5);
    c = (c-1)/(alpha); % -1..1 in the transition region
    c=min(c, 1);
    c=max(c, -1);
    RRC_h = sqrt(1/2+cos(pi/2*(c+1))/2);
    
    ir = real(ifft(RRC_h));    
end

% ************************************
% remove any energy at frequencies > BW/2
% ************************************
function s = brickwallFilter(s, rate, BW)
    n = numel(s);
    fb = FFT_frequencyBasis(n, rate);
    ix = find(abs(fb) > BW / 2);
    s = fft(s);
    s(ix) = 0;
    s = real(ifft(s));
end

% ************************************
% for an impulse response s at rate, write the
% frequency response to fname
% ************************************
function exportFrequencyResponse(s, rate, fname)
    fb = fftshift(FFT_frequencyBasis(numel(s), rate));
    fr = fftshift(fft(s));
    figure(335); grid on;
    plot(fb, 20*log10(abs(fr)));
    title('exported frequency response');
    xlabel('normalized frequency');
    ylabel('dB');
    save(fname, 'fb', 'fr');
end

% ************************************
% calculates the frequency that corresponds to
% each FFT bin (negative, zero, positive)
% ************************************
function fb_Hz = FFT_frequencyBasis(n, rate_Hz)
    fb = 0:(n - 1);
    fb = fb + floor(n / 2);
    fb = mod(fb, n);
    fb = fb - floor(n / 2);
    fb = fb / n; % now [0..0.5[, [-0.5..0[
    fb_Hz = fb * rate_Hz;
end

% *******************************************************
% delay-matching between two signals (complex/real-valued)
%
% * matches the continuous-time equivalent waveforms
%   of the signal vectors (reconstruction at Nyquist limit =>
%   ideal lowpass filter)
% * Signals are considered cyclic. Use arbitrary-length 
%   zero-padding to turn a one-shot signal into a cyclic one.
%
% * output:
%   => coeff: complex scaling factor that scales 'ref' into 'signal'
%   => delay 'deltaN' in units of samples (subsample resolution)
%      apply both to minimize the least-square residual      
%   => 'shiftedRef': a shifted and scaled version of 'ref' that 
%      matches 'signal' 
%   => (signal - shiftedRef) gives the residual (vector error)
%
% Example application
% - with a full-duplex soundcard, transmit an arbitrary cyclic test signal 'ref'
% - record 'signal' at the same time
% - extract one arbitrary cycle
% - run fitSignal
% - deltaN gives the delay between both with subsample precision
% - 'shiftedRef' is the reference signal fractionally resampled 
%   and scaled to optimally match 'signal'
% - to resample 'signal' instead, exchange the input arguments
% *******************************************************
function [coeff, shiftedRef, deltaN] = fitSignal_FFT(signal, ref)
    n=length(signal);
    % xyz_FD: Frequency Domain
    % xyz_TD: Time Domain
    % all references to 'time' and 'frequency' are for illustration only

    forceReal = isreal(signal) && isreal(ref);
    
    % *******************************************************
    % Calculate the frequency that corresponds to each FFT bin
    % [-0.5..0.5[
    % *******************************************************
    binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;

    % *******************************************************
    % Delay calculation starts:
    % Convert to frequency domain...
    % *******************************************************
    sig_FD = fft(signal);
    ref_FD = fft(ref, n);

    % *******************************************************
    % ... calculate crosscorrelation between 
    % signal and reference...
    % *******************************************************
    u=sig_FD .* conj(ref_FD);
    if mod(n, 2) == 0
        % for an even sized FFT the center bin represents a signal
        % [-1 1 -1 1 ...] (subject to interpretation). It cannot be delayed. 
        % The frequency component is therefore excluded from the calculation.
        u(length(u)/2+1)=0;
    end
    Xcor=abs(ifft(u));

    %  figure(); plot(abs(Xcor));
    
    % *******************************************************
    % Each bin in Xcor corresponds to a given delay in samples.
    % The bin with the highest absolute value corresponds to
    % the delay where maximum correlation occurs.
    % *******************************************************
    integerDelay = find(Xcor==max(Xcor));
    
    % (1): in case there are several bitwise identical peaks, use the first one
    % Minus one: Delay 0 appears in bin 1
    integerDelay=integerDelay(1)-1;

    % Fourier transform of a pulse shifted by one sample
    rotN = exp(2i*pi*integerDelay .* binFreq);

    uDelayPhase = -2*pi*binFreq;
    
    % *******************************************************
    % Since the signal was multiplied with the conjugate of the
    % reference, the phase is rotated back to 0 degrees in case
    % of no delay. Delay appears as linear increase in phase, but
    % it has discontinuities.
    % Use the known phase (with +/- 1/2 sample accuracy) to 
    % rotate back the phase. This removes the discontinuities.
    % *******************************************************
    %  figure(); plot(angle(u)); title('phase before rotation');
    u=u .* rotN;
    
    % figure(); plot(angle(u)); title('phase after rotation');
    
    % *******************************************************
    % Obtain the delay using linear least mean squares fit
    % The phase is weighted according to the amplitude.
    % This suppresses the error caused by frequencies with
    % little power, that may have radically different phase.
    % *******************************************************
    weight = abs(u); 
    constRotPhase = 1 .* weight;
    uDelayPhase = uDelayPhase .* weight;
    ang = angle(u) .* weight;
    r = [constRotPhase; uDelayPhase] .' \ ang.'; %linear mean square
    
    %rotPhase=r(1); % constant phase rotation, not used.
    % the same will be obtained via the phase of 'coeff' further down
    fractionalDelay=r(2);
    
    % *******************************************************
    % Finally, the total delay is the sum of integer part and
    % fractional part.
    % *******************************************************
    deltaN = integerDelay + fractionalDelay;

    % *******************************************************
    % provide shifted and scaled 'ref' signal
    % *******************************************************
    % this is effectively time-convolution with a unit pulse shifted by deltaN
    rotN = exp(-2i*pi*deltaN .* binFreq);
    ref_FD = ref_FD .* rotN;
    shiftedRef = ifft(ref_FD);
    
    % *******************************************************
    % Again, crosscorrelation with the now time-aligned signal
    % *******************************************************
    coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
    shiftedRef=shiftedRef * coeff;

    if forceReal
        shiftedRef = real(shiftedRef);
    end
end

Baseband-equivalent phase noise model

Markus Nentwig December 18, 20114 comments Coded in Matlab
% ****************************************************************
% baseband equivalent source of local oscillator with phase noise
% Markus Nentwig, 18.12.2011
% ****************************************************************
function pn_generator()
    close all;
    
    % ****************************************************************
    % PN generator configuration
    % ****************************************************************
    srcPar = struct();
    srcPar.n = 2 ^ 18; % generated number of output samples
    srcPar.rate_Hz = 7.68e6; % sampling rate
    srcPar.f_Hz = [0, 10e3, 1e6, 9e9]; % phase noise spectrum, frequencies
    srcPar.g_dBc1Hz = [-80, -80, -140, -140]; % phase noise spectrum, magnitude
    srcPar.spursF_Hz = [300e3, 400e3, 700e3]; % discrete spurs (set [] if not needed)
    srcPar.spursG_dBc = [-50, -55, -60]; % discrete spurs, power relative to carrier
    
    % ****************************************************************
    % run PN generator
    % ****************************************************************
    s = PN_src(srcPar);
    
    if false
        % ****************************************************************
        % export phase noise baseband-equivalent signal for use in simulator etc
        % ****************************************************************
        tmp = [real(s); imag(s)] .';
        save('phaseNoiseSample.dat', 'tmp', '-ascii');
    end
    
    if exist('spectrumAnalyzer', 'file')
        % ****************************************************************
        % spectrum analyzer configuration
        % ****************************************************************
        SAparams = struct();
        SAparams.rate_Hz = srcPar.rate_Hz; % sampling rate of the input signal
        SAparams.pRef_W = 1e-3; % unity signal represents 0 dBm (1/1000 W)
        SAparams.pNom_dBm = 0; % show 0 dBm as 0 dB;
        SAparams.filter = 'brickwall';
        SAparams.RBW_window_Hz = 1000; % convolve power spectrum with a 1k filter
        SAparams.RBW_power_Hz = 1; % show power density as dBc in 1 Hz
        SAparams.noisefloor_dB = -250; % don't add artificial noise
        SAparams.logscale = true; % use logarithmic frequency axis
        
        % plot nominal spectrum
        figure(1); grid on; 
        h = semilogx(max(srcPar.f_Hz, 100) / 1e6, srcPar.g_dBc1Hz, 'k+-');
        set(h, 'lineWidth', 3);
        hold on;
        spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'fig', 1);
        ylabel('dBc in 1 Hz');
        legend('nominal PSD', 'output spectrum');
        title('check match with nominal PSD');
        spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'RBW_power_Hz', 'sine', 'fig', 2);
        title('check carrier level (0 dBc); check spurs level(-50/-55/-60 dBc)');
        ylabel('dBc for continuous-wave signal');
        spectrumAnalyzer('signal', s, SAparams, 'fig', 3, 'logscale', false);
        ylabel('dBc in 1 Hz');    
    end
end

function pn_td = PN_src(varargin)
    def = {'includeCarrier', true, ...
           'spursF_Hz', [], ...
           'spursG_dBc', [], ...
           'fMax_Hz', []};
    p = vararginToStruct(def, varargin);
    
    % length of signal in the time domain (after ifft)
    len_s = p.n / p.rate_Hz

    % FFT bin frequency spacing
    deltaF_Hz = 1 / len_s
    
    % construct AWGN signal in the frequency domain     
    % a frequency domain bin value of n gives a time domain power of 1
    % for example ifft([4 0 0 0]) => 1 1 1 1 
    % each bin covers a frequency interval of deltaF_Hz
    mag = p.n;
    
    % scale "unity power in one bin" => "unity power per Hz":
    % multiply with sqrt(deltaF_Hz): 
    mag = mag * sqrt(deltaF_Hz);
    
    % Create noise according to mag in BOTH real- and imaginary value
    mag = mag * sqrt(2);

    % both real- and imaginary part contribute unity power => divide by sqrt(2)
    pn_fd = mag / sqrt(2) * (randn(1, p.n) + 1i * randn(1, p.n));
    
    % frequency vector corresponding to the FFT bins (0, positive, negative)
    fb_Hz = FFT_freqbase(p.n, deltaF_Hz);
    
    % interpolate phase noise spectrum on frequency vector
    % note: interpolate dB on logarithmic frequency axis
    H_dB = interp1(log(p.f_Hz+eps), p.g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');

    % dB => magnitude
    H = 10 .^ (H_dB / 20);
    %    H = 1e-6; % sanity check: enforce flat -120 dBc in 1 Hz
    
    % apply filter to noise spectrum
    pn_fd = pn_fd .* H;

    % set spurs
    for ix = 1:numel(p.spursF_Hz)
        fs = p.spursF_Hz(ix);
        u = abs(fb_Hz - fs);
        ix2 = find(u == min(u), 1);

        % random phase
        rot = exp(2i*pi*rand());
        
        % bin value of n: unity (carrier) power (0 dBc)
        % scale with sqrt(2) because imaginary part will be discarded
        % scale with sqrt(2) because the tone appears at both positive and negative frequencies
        smag = 2 * p.n * 10 ^ (p.spursG_dBc(ix) / 20);
        pn_fd(ix2) = smag * rot;
    end
    
    % limit bandwidth (tool to avoid aliasing in an application
    % using the generated phase noise signal)
    if ~isempty(p.fMax_Hz)
        pn_fd(find(abs(fb_Hz) > p.fMax_Hz)) = 0;
    end
    
    % convert to time domain
    pn_td = ifft(pn_fd);

    % discard imaginary part 
    pn_td = real(pn_td);

    % Now pn_td is a real-valued random signal with a power spectral density 
    % as specified in f_Hz / g_dBc1Hz.
    
    % phase-modulate to carrier
    % note: d/dx exp(x) = 1 near |x| = 1
    % in other words, the phase modulation has unity gain for small phase error
    pn_td = exp(i*pn_td);
    
    if ~p.includeCarrier
        % remove carrier
        % returns isolated phase noise component
        pn_td = pn_td - 1;
    end
end    

% returns a vector of frequencies corresponding to n FFT bins, when the
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
    fb_Hz = 0:(n - 1);
    fb_Hz = fb_Hz + floor(n / 2);
    fb_Hz = mod(fb_Hz, n);
    fb_Hz = fb_Hz - floor(n / 2);
    fb_Hz = fb_Hz * deltaF_Hz;
end

% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
%   varargin = {'hello'}
    r = flattenCellArray(varargin, struct());
end

function r = flattenCellArray(arr, r)
    ix=1;
    ixMax = numel(arr);
    while ix <= ixMax
        e = arr{ix};
        
        if iscell(e)
            % cell array at 'key' position gets recursively flattened
            % becomes struct
            r = flattenCellArray(e, r);
        elseif ischar(e)
            % string => key.
            % The following entry is a value
            ix = ix + 1;
            v = arr{ix};
            % store key-value pair
            r.(e) = v;
        elseif isstruct(e)
            names = fieldnames(e);
            for ix2 = 1:numel(names)
                k = names{ix2};
                r.(k) = e.(k);
            end
        else
            e
            assert(false)
        end
        ix=ix+1;
    end % while
end

Digital model for analog filters

Markus Nentwig December 17, 2011 Coded in Matlab
% discrete-time model for Laplace-domain expression
% Markus Nentwig, 30.12.2011
function sn_model()
    close all;
    run_demo1(10);
    run_demo2(11);
end

% ************************************
% Constructs a FIR model for a relatively
% narrow-band continuous-time IIR filter.
% At the edge of the first Nyquist zone 
% (+/- fSample/2), the frequency response
% is down by about 70 dB, which makes 
% the modeling unproblematic. 
% The impact of different windowing options 
% is visible both at high frequencies, but
% also as deviation between original frequency
% response and model at the passband edge.
% ************************************
function run_demo1(fig)
    [b, a] = getContTimeExampleFilter();
    fc_Hz = 0.5e6; % frequency corresponding to omegaNorm == 1    
    
    commonParameters = struct(...
        's_a', a, ...
        's_b', b, ...
        'z_rate_Hz', 3e6, ...
        's_fNorm_Hz', fc_Hz, ...
        'fig', fig);
    
    % sample impulse response without windowing
    ir = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'plotstyle_s', 'k-', ...
        'plotstyle_z', 'b-');

    % use mild windowing
    ir = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'plotstyle_z', 'r-', ...
        'winLen_percent', 4);

    % use heavy windowing - error shows at passband edge
    ir = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'plotstyle_z', 'm-', ...
        'winLen_percent', 100);
    
    legend('s-domain', 'sampled, no window', 'sampled, 4% RC window', 'sampled, 100% RC window')
    
    figure(33); clf;
    h = stem(ir);
    set(h, 'markersize', 2);
    set(h, 'lineWidth', 2);
    title('sampled impulse response');
end

% ************************************
% model for a continuous-time IIR filter
% that is relatively wide-band.
% The frequency response requires some 
% manipulation at the edge of the Nyquist zone.
% Otherwise, there would be an abrupt change
% that would result in an excessively long impulse
% response.
% ************************************
function run_demo2(fig)
    [b, a] = getContTimeExampleFilter();
    fc_Hz = 1.4e6; % frequency corresponding to omegaNorm == 1    
    
    commonParameters = struct(...
        's_a', a, ...
        's_b', b, ...
        'z_rate_Hz', 3e6, ...
        's_fNorm_Hz', fc_Hz, ...
        'fig', fig);
    
    % sample impulse response without any manipulations
    ir1 = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'aliasZones', 0, ...
        'plotstyle_s', 'k-', ...
        'plotstyle_z', 'b-');

    % use artificial aliasing (introduces some passband error)
    ir2 = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'aliasZones', 5, ...
        'plotstyle_z', 'r-');

    % use artificial low-pass filter (freq. response
    % becomes invalid beyond +/- BW_Hz / 2)
    ir3 = sampleLaplaceDomainImpulseResponse(...
        commonParameters, ...
        'aliasZones', 0, ...
        'BW_Hz', 2.7e6, ...
        'plotstyle_z', 'm-');
    line([0, 2.7e6/2, 2.7e6/2], [-10, -10, -50]);
    
    legend('s-domain', ...
           sprintf('unmodified (%i taps)', numel(ir1)), ...
           sprintf('artificial aliasing (%i taps)', numel(ir2)), ...
           sprintf('artificial LP filter (%i taps)', numel(ir3)));
    title('2nd example: wide-band filter model frequency response');
    
    figure(350); grid on; hold on;
    subplot(3, 1, 1);
    stem(ir1, 'b'); xlim([1, numel(ir1)])
    title('wide-band filter model: unmodified');
    subplot(3, 1, 2);
    stem(ir2, 'r');xlim([1, numel(ir1)]);
    title('wide-band filter model: art. aliasing');
    subplot(3, 1, 3);
    stem(ir3, 'm');xlim([1, numel(ir1)]);
    title('wide-band filter model: art. LP filter');
end

% Build example Laplace-domain model
function [b, a] = getContTimeExampleFilter()
    if true
        order = 6;
        ripple_dB = 1.2;
        omegaNorm = 1;
        [b, a] = cheby1(order, ripple_dB, omegaNorm, 's');
    else
        % same as above, if cheby1 is not available
        b =  0.055394;
        a = [1.000000   0.868142   1.876836   1.109439   0.889395   0.279242   0.063601];
    end
end

% ************************************
% * Samples the impulse response of a Laplace-domain
%   component
% * Adjusts group delay and impulse response length so that
%   the discarded part of the impulse response is
%   below a threshold.
% * Applies windowing 
% 
% Mandatory arguments:
% s_a, s_b:    Laplace-domain coefficients
% s_fNorm_Hz:  normalization frequency for a, b coefficients
% z_rate_Hz:   Sampling rate of impulse response
% 
% optional arguments:
% truncErr_dB: Maximum truncation error of impulse response
% nPts:        Computed impulse response length before truncation
% winLen_percent: Part of the IR where windowing is applied
% BW_Hz:       Apply artificical lowpass filter for |f| > +/- BW_Hz / 2
% 
% plotting:
% fig:         Figure number
% plotstyle_s: set to plot Laplace-domain frequency response
% plotstyle_z: set to plot z-domain model freq. response
% ************************************
function ir = sampleLaplaceDomainImpulseResponse(varargin)
    def = {'nPts', 2^18, ...
           'truncErr_dB', -60, ...
           'winLen_percent', -1, ...
           'fig', 99, ...
           'plotstyle_s', [], ...
           'plotstyle_z', [], ...
           'aliasZones', 1, ...
           'BW_Hz', []};
    p = vararginToStruct(def, varargin);

    % FFT bin frequencies
    fbase_Hz = FFT_frequencyBasis(p.nPts, p.z_rate_Hz);
    
    % instead of truncating the frequency response at +/- z_rate_Hz, 
    % fold the aliases back into the fundamental Nyquist zone.
    % Otherwise, we'd try to build a near-ideal lowpass filter,
    % which is essentially non-causal and requires a long impulse
    % response with artificially introduced group delay.
    H = 0;
    for alias = -p.aliasZones:p.aliasZones
        
        % Laplace-domain "s" in normalized frequency
        s = 1i * (fbase_Hz + alias * p.z_rate_Hz) / p.s_fNorm_Hz;
        
        % evaluate polynomial in s
        H_num = polyval(p.s_b, s);
        H_denom = polyval(p.s_a, s);
        Ha = H_num ./ H_denom;
        H = H + Ha;
    end
    
    % plot |H(f)| if requested
    if ~isempty(p.plotstyle_s)
        figure(p.fig); hold on; grid on;
        fr = fftshift(20*log10(abs(H) + eps));
        h = plot(fftshift(fbase_Hz), fr, p.plotstyle_s);
        set(h, 'lineWidth', 3);
        xlabel('f/Hz'); ylabel('|H(f)| / dB');
        xlim([0, p.z_rate_Hz/2]);
    end

    % apply artificial lowpass filter
    if ~isempty(p.BW_Hz)
        % calculate RC transition bandwidth
        BW_RC_trans_Hz = p.z_rate_Hz - p.BW_Hz;
        
        % alpha (RC filter parameter) implements the 
        % RC transition bandwidth:
        alpha_RC = BW_RC_trans_Hz / p.z_rate_Hz / 2;

        % With a cutoff frequency of BW_RC, the RC filter 
        % reaches |H(f)| = 0 at f = z_rate_Hz / 2
        % BW * (1+alpha) = z_rate_Hz / 2
        BW_RC_Hz = p.z_rate_Hz / ((1+alpha_RC));
        HRC = raisedCosine(fbase_Hz, BW_RC_Hz, alpha_RC);
        H = H .* HRC;
    end
        
    % frequency response => impulse response
    ir = ifft(H);
    
    % assume s_a, s_b describe a real-valued impulse response
    ir = real(ir);
    
    % the impulse response peak is near the first bin.
    % there is some earlier ringing, because evaluating the s-domain
    % expression only for s < +/- z_rate_Hz / 2 implies an ideal, 
    % non-causal low-pass filter.

    ir = fftshift(ir);
    % now the peak is near the center
    
    % threshold for discarding samples
    peak = max(abs(ir));    
    thr = peak * 10 ^ (p.truncErr_dB / 20);

    % first sample that is above threshold
    % determines the group delay of the model
    ixFirst = find(abs(ir) > thr, 1, 'first');

    % last sample above threshold
    % determines the length of the impulse response
    ixLast = find(abs(ir) > thr, 1, 'last');
    
    % truncate
    ir = ir(ixFirst:ixLast);

    % apply windowing
    if p.winLen_percent > 0
        % note: The window drops to zero for the first sample that
        % is NOT included in the impulse response.
        v = linspace(-1, 0, numel(ir)+1);
        v = v(1:end-1);
        v = v * 100 / p.winLen_percent;
        v = v + 1;
        v = max(v, 0);
        win = (cos(v*pi) + 1) / 2;
        ir = ir .* win;
    end
    
    % plot frequency response, if requested
    if ~isempty(p.plotstyle_z)
        irPlot = zeros(1, p.nPts);
        irPlot(1:numel(ir)) = ir;
        figure(p.fig); hold on;
        fr = fftshift(20*log10(abs(fft(irPlot)) + eps));
        h = plot(fftshift(fbase_Hz), fr, p.plotstyle_z);
        set(h, 'lineWidth', 3);
        xlabel('f/Hz'); ylabel('|H(f)| / dB');
        xlim([0, p.z_rate_Hz/2]);
    end
end

% ************************************
% raised cosine frequency response
% ************************************
function H = raisedCosine(f, BW, alpha)
    c=abs(f * 2 / BW);
    
    % c=-1 for lower end of transition region
    % c=1 for upper end of transition region
    c=(c-1) / alpha;
    
    % clip to -1..1 range
    c=min(c, 1);
    c=max(c, -1);
    
    H = 1/2+cos(pi/2*(c+1))/2;
end

% ************************************
% calculates the frequency that corresponds to each
% FFT bin (negative, zero, positive)
% ************************************
function fb_Hz = FFT_frequencyBasis(n, rate_Hz)
    fb = 0:(n - 1);
    fb = fb + floor(n / 2);
    fb = mod(fb, n);
    fb = fb - floor(n / 2);
    fb = fb / n; % now [0..0.5[, [-0.5..0[
    fb_Hz = fb * rate_Hz;
end

% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
%   varargin = {'hello'}
    r = flattenCellArray(varargin, struct());
end

function r = flattenCellArray(arr, r)
    ix=1;
    ixMax = numel(arr);
    while ix <= ixMax
        e = arr{ix};
        
        if iscell(e)
            % cell array at 'key' position gets recursively flattened
            % becomes struct
            r = flattenCellArray(e, r);
        elseif ischar(e)
            % string => key.
            % The following entry is a value
            ix = ix + 1;
            v = arr{ix};
            % store key-value pair
            r.(e) = v;
        elseif isstruct(e)
            names = fieldnames(e);
            for ix2 = 1:numel(names)
                k = names{ix2};
                r.(k) = e.(k);
            end
        else
            e
            assert(false)
        end
        ix=ix+1;
    end % while
end

Testing the Flat-Top Windowing Function

Rick Lyons December 14, 2011 Coded in Matlab
%    Code for testing the 'Wind_Flattop(Spec)' function in 
%    reducing 'scalloping loss' errors in time signal amplitude 
%    estimation.
%
%    Generates a time-domain sinusoid, computes its FFT, and passes 
%    that FFT sequence to the 'Wind_Flattop(Spec)' function.
%
%    The maximum output sample of the 'Wind_Flattop(Spec)' function
%    is used to estimate the peak amplitude of the original 
%    sinusoidal time-domain test signal.
%
%    The user controls the values for the test sinusoid's  
%    'Test_Freq' and peak amplitude 'Peak_Amp' near lines 23 & 24. 
%
%    Richard Lyons [December, 2011]

clear all, clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define test parameters
Test_Freq = 7.22; % Test tone's freq. Must be less that N/2
Peak_Amp = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 64; % Number of time samples
Index = (0:N-1);

X = Peak_Amp*cos(2*pi*(Test_Freq)*Index/N + pi/3);

figure(1), clf
subplot(2,1,1)
plot(X,':ko', 'markersize', 4)
title('Original time signal'), grid on, zoom on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   FFT the input 'X' sequence and call 'Wind_Flattop()' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Spec = fft(X);
[Windowed_Spec] = Wind_Flattop(Spec);

subplot(2,1,2), plot(abs(Spec),'ko', 'markersize', 3)
title('SpecMag of unwindowed time signal'), grid on, zoom on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Display results accuracy (error)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(['Test Freq = ',num2str(Test_Freq),...
        ',  True Peak Amplitude = ',num2str(Peak_Amp)])

Mag_peak_unwindowed = max(abs(Spec));
Unwindowed_Amp_Estimate = 2*Mag_peak_unwindowed/N;
Unwindowed_Amp_Estimate_Error_in_dB = ...
    20*log10(Unwindowed_Amp_Estimate/Peak_Amp);
disp(' ')
disp(['Unwindowed Peak Amplitude Estimate = ',...
        num2str(Unwindowed_Amp_Estimate)])
disp(['Unwindowed Estimate Error in dB = ',...
        num2str(Unwindowed_Amp_Estimate_Error_in_dB),' dB'])

M_peak_windowed = max(abs(Windowed_Spec));
Windowed_Amp_Estimated = 2*M_peak_windowed/N;
Windowed_Amp_Estimation_Error_in_dB = ...
    20*log10(Windowed_Amp_Estimated/Peak_Amp);
disp(' ')
disp(['Windowed Peak Amplitude Estimate = ',...
        num2str(Windowed_Amp_Estimated)])
disp(['Windowed Estimate Error in dB = ',...
        num2str(Windowed_Amp_Estimation_Error_in_dB),' dB'])

Flat-Top Windowing Function for the Accurate Measurement of a Sinusoid's Peak Amplitude Based on FFT Data

Rick Lyons December 14, 2011 Coded in Matlab
function [Windowed_Spec] = Wind_Flattop(Spec)

%  Given an input spectral sequence 'Spec', that is the 
%  FFT of some time sequence 'x', Wind_Flattop(Spec) 
%  returns a spectral sequence that is equivalent
%  to the FFT of a flat-top windowed version of time 
%  sequence 'x'.  The peak magnitude values of output 
%  sequence 'Windowed_Spec' can be used to accurately 
%  estimate the peak amplitudes of sinusoidal components 
%  in time sequence 'x'.

%   Input: 'Spec' (a sequence of complex FFT sample values)
%
%   Output: 'Windowed_Spec' (a sequence of complex flat-top  
%                            windowed FFT sample values)
%
%   Based on Lyons': "Reducing FFT Scalloping Loss Errors 
%   Without Multiplication", IEEE Signal Processing Magazine, 
%   DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
%   Richard Lyons [December, 2011]

N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];

% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
    +g_Coeffs(2)*Spec(N) + Spec(1) ...
    +g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);

Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
    +g_Coeffs(2)*Spec(1) + Spec(2) ...
    +g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);

% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
    +g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
    +g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);

Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
    +g_Coeffs(2)*Spec(N-1) + Spec(N) ...
    +g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);

% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
	Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
        +g_Coeffs(2)*Spec(K-1) + Spec(K) ...
        +g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function

16QAM Modem model (Basic)

Kadhiem Ayob December 12, 2011 Coded in Matlab
clear all; close all;

n = 2^16;   %number of symbols
Fsym = 12.5; %Msps
Fs = 100;   %MHz, IF sampling frequency
Fc = 20;    %MHz, upconverter frequency

%generate 16QAM raw symbols
alphabet = [-1 -1/3 +1/3 +1];
tx = complex(randsrc(1,n,alphabet),randsrc(1,n,alphabet));

%pulse shaping, root raised cosine
h = firrcos(50,Fsym/4,.15,Fsym,'rolloff','sqrt');
tx_up1 = zeros(1,2*n);
tx_up1(1:2:end-1) = tx;
tx_shaped = filter(2*h,1,tx_up1);

%further upsampling by 4
tx_up2 = resample(tx_shaped,4,1);

%upconvert to 10MHz centre
f1 = exp(1i*2*pi*(0:8*n-1)*Fc/Fs);
tx_upconverted = tx_up2 .* f1;

%channel signal
rx_real = real(tx_upconverted);

%Rx shifts signal back to zero 
f2 = exp(1i*2*pi*(0:8*n-1)*-Fc/Fs);
rx_downconverted = rx_real.*f2;

%Rx downsamples back by 4
rx_dn1 = resample(rx_downconverted,1,4);

%rx matched filter
rx_dn2 = filter(2*h,1,rx_dn1);
rx_dn2 = rx_dn2(50:end); %remove initial zeros

%one phase should be correct(odd/even)
rx_e = rx_dn2(1:2:end); %odd
rx_o = rx_dn2(2:2:end); %even

subplot(3,1,1);plot(real(tx),imag(tx),'o');
axis([-2 +2 -2 +2]);grid;
title('Tx constellations')
subplot(3,1,2);plot(real(rx_e),imag(rx_e),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, first phase')
subplot(3,1,3);plot(real(rx_o),imag(rx_o),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, second phase')

Phase drift of NCO

Kadhiem Ayob December 10, 2011 Coded in Matlab
clear all; close all;

n = 50000;                            %number of test samples
Fs = 100;                             %sampling rate in Msps
Fo = 1.17;                            %target NCO frequency in MHz
M = [2^16 10000];                     %NCO phase resolution, two cases
for i = 1:2
    inc = round(M(i)*Fo/Fs);                         %phase increment
    lut = round(2047*cos(2*pi*(0:M(i)-1)/M(i)));     %lut, one cycle cosine data
    addr = 0:inc:inc*(n-1);
    addr = round(mod(addr,M(i)));
    addr(addr >= M(i)) = addr(addr >= M(i)) - M(i);  %check address overflow
    y(i,1:n) = lut(addr+1);                          %add 1 for matlab LUT
end

fprintf('increment value = %2.6f\r',M*Fo/Fs);

subplot(2,1,1);hold;
plot([y(1,1:1000) zeros(1,100) y(1,end-999:end)],'.-');
plot([y(2,1:1000) zeros(1,100) y(2,end-999:end)],'r.--');
legend('drift','no drift');
title('initial and last 1000 samples');

[P1, F1] = pwelch(y(1,:), hann(n), 0, n, Fs);
[P2, F2] = pwelch(y(2,:), hann(n), 0, n, Fs);
subplot(2,1,2); hold
plot(F1,10*log10(P1));
plot(F2,10*log10(P2),'r--');
legend('drift','no drift')
axis ([0 50 -20 100]);
grid
xlabel('MHz')

Spectrum analyzer

Markus Nentwig December 10, 20111 comment Coded in Matlab
% ************************************
% spectrum analyzer
% Markus Nentwig, 10.12.2011
%
% Note: If the signal is not cyclic, apply conventional windowing before 
% calling spectrumAnalyzer
% ************************************

function satest()
    close all;
    rate_Hz = 30.72e6; % sampling rate
    n = 100000; % number of samples

    noise = randn(1, n);
    pulse = zeros(1, n); pulse(1) = 1;
    
    % ************************************
    % example (linear frequency scale, RBW filter)
    % ************************************
    tmp1 = normalize(RRC_filter('signal', noise, 'rate_Hz', rate_Hz, 'alpha', 0.22, 'BW_Hz', 3.84e6, 'fCenter_Hz', 5e6));

    % pRef = 0.001: a normalized signal will show as 0 dBm = 0.001 W
    saPar = struct('rate_Hz', rate_Hz, 'fig', 1, 'pRef_W', 1e-3, 'RBW_power_Hz', 3.84e6, 'filter', 'brickwall', 'plotStyle', 'b-'); 
    spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 1e3);
    spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 30e3, 'plotStyle', 'k-');
    spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 300e3, 'plotStyle', 'r-');
    spectrumAnalyzer(saPar, 'signal', tmp1, 'RBW_window_Hz', 30e3, 'filter', 'gaussian', 'plotStyle', 'g-');
    legend('1k RBW filter', '30k RBW filter', '300k RBW filter', '30k Gaussian filter (meas. instrument model)');
    title('WCDMA-like signal');
    
    % ************************************
    % example (logarithmic frequency scale, 
    % no RBW filter)
    % ************************************
    fPar = struct('rate_Hz', rate_Hz, ...
                  'chebyOrder', 6, ...
                  'chebyRipple_dB', 1, ...
                  'fCorner_Hz', 1e5);

    saPar = struct('rate_Hz', rate_Hz, ...
                   'pRef_W', 1e-3, ...
                   'fMin_Hz', 1e3, ...
                   'RBW_power_Hz', 1e5, ...
                   'RBW_window_Hz', 1e3, ...
                   'filter', 'none', ...
                   'plotStyle', 'b-', ...
                   'logscale', true, ...
                   'fig', 2, ...
                   'noisefloor_dB', -300);
    
    tmp1 = normalize(IIR_filterExample('signal', noise, fPar));
    tmp2 = normalize(IIR_filterExample('signal', pulse, fPar));
    
    spectrumAnalyzer('signal', tmp1, saPar);
    spectrumAnalyzer('signal', tmp2, saPar, 'plotStyle', 'k-');
end

function sig = normalize(sig)
    sig = sig / sqrt(sig * sig' / numel(sig));
end

% ************************************
% calculates the frequency that corresponds to each
% FFT bin (negative, zero, positive)
% ************************************
function fb_Hz = FFT_frequencyBasis(n, rate_Hz)
    fb = 0:(n - 1);
    fb = fb + floor(n / 2);
    fb = mod(fb, n);
    fb = fb - floor(n / 2);
    fb = fb / n; % now [0..0.5[, [-0.5..0[
    fb_Hz = fb * rate_Hz;
end

% ************************************
% root-raised cosine filter (to generate
% example signal)
% ************************************
function sig = RRC_filter(varargin)
    def = struct('fCenter_Hz', 0);
    p = vararginToStruct(def, varargin);
    
    n = numel(p.signal);
    fb_Hz = FFT_frequencyBasis(n, p.rate_Hz);

    % frequency relative to cutoff frequency
    c = abs((fb_Hz - p.fCenter_Hz) / (p.BW_Hz / 2));
    % c=-1 for lower end of transition region
    % c=1 for upper end of transition region
    c=(c-1) / p.alpha;

    % clip to -1..1 range
    c=min(c, 1);
    c=max(c, -1);
    
    % raised-cosine filter
    H = 1/2+cos(pi/2*(c+1))/2;

    % root-raised cosine filter
    H = sqrt(H);

    % apply filter
    sig = ifft(fft(p.signal) .* H);
end    

% ************************************
% continuous-time filter example
% ************************************
function sig = IIR_filterExample(varargin)
    p = vararginToStruct(varargin);
    
    % design continuous-time IIR filter
    [IIR_b, IIR_a] = cheby1(p.chebyOrder, p.chebyRipple_dB, 1, 's');
    
    % evaluate on the FFT frequency grid
    fb_Hz = FFT_frequencyBasis(numel(p.signal), p.rate_Hz);
    
    % Laplace domain operator, normalized to filter cutoff frequency
    % (the above cheby1 call designs for unity cutoff)
    s = 1i * fb_Hz / p.fCorner_Hz;
    
    % polynomial in s
    H_num = polyval(IIR_b, s);
    H_denom = polyval(IIR_a, s);
    H = H_num ./ H_denom;
    
    % apply filter
    sig = ifft(fft(p.signal) .* H);
end

% ----------------------------------------------------------------------
% optionally: Move the code below into spectrumAnalyzer.m
function [fr, fb, handle] = spectrumAnalyzer(varargin)
    def = struct(...
        'noisefloor_dB', -150, ...
        'filter', 'none', ...
        'logscale', false,  ...
        'NMax', 10000, ...
        'freqUnit', 'MHz', ...
        'fig', -1, ...
        'plotStyle', 'b-', ...
        'originMapsTo_Hz', 0 ...
        );

    p = vararginToStruct(def, varargin);
    signal = p.signal;

    handle=nan; % avoid warning
    
    % A resolution bandwidth value of 'sine' sets the RBW to the FFT bin spacing.
    % The power of a pure sine wave now shows correctly from the peak in the spectrum (unity => 0 dB)
    singleBinMode=strcmp(p.RBW_power_Hz, 'sine');

    nSamples = numel(p.signal);
    binspacing_Hz = p.rate_Hz / nSamples;
    windowBW=p.RBW_window_Hz;
    noisefloor=10^(p.noisefloor_dB/10);
    % factor in the scaling factor that will be applied later for conversion to
    % power in RBW
    if ~singleBinMode
        noisefloor = noisefloor * binspacing_Hz / p.RBW_power_Hz;
    end

    % fr: "f"requency "r"esponse (plot y data)
    % fb: "f"requency "b"ase (plot x data)
    fr = fft(p.signal);
    
    scale_to_dBm=sqrt(p.pRef_W/0.001);

    % Normalize amplitude to the number of samples that contribute
    % to the spectrum
    fr=fr*scale_to_dBm/nSamples;

    % magnitude
    fr = fr .* conj(fr);

    [winLeft, winRight] = spectrumAnalyzerGetWindow(p.filter, singleBinMode, p.RBW_window_Hz, binspacing_Hz);
    % winLeft is the right half of the window, but it appears on the
    % left side of the FFT space

    winLen=0;
    if ~isempty(winLeft)  
        
        % zero-pad the power spectrum in the middle with a length
        % equivalent to the window size. 
        % this guarantees that there is enough bandwidth for the filter!
        % one window length is enough, the spillover from both sides overlaps
        % Scale accordingly.
        winLen=size(winLeft, 2)+size(winRight, 2);
        
        % note: fr is the power shown in the plot, NOT a frequency
        % domain representation of a signal.
        % There is no need to renormalize because of the length change
        center=floor(nSamples/2)+1;
        rem=nSamples-center;
        fr=[fr(1:center-1), zeros(1, winLen-1), fr(center:end)];
        % construct window with same length as fr
        win=zeros(size(fr));
        win(1:1+size(winLeft, 2)-1)=winLeft;
        win(end-size(winRight, 2)+1:end)=winRight;
        
        assert(isreal(win));
        assert(isreal(fr));
        assert(size(win, 2)==size(fr, 2));
        
        % convolve using FFT
        win=fft(win);
        fr=fft(fr);
        
        fr=fr .* win;
        fr=ifft(fr);
        fr=real(fr); % chop off roundoff error imaginary part
        fr=max(fr, 0); % chop off small negative numbers
        
        % remove padding
        fr=[fr(1:center-1), fr(end-rem+1:end)];
    end

    % ************************************
    % build frequency basis and rotate 0 Hz to center
    % ************************************
    fb = FFT_frequencyBasis(numel(fr), p.rate_Hz);    
    fr = fftshift(fr);
    fb = fftshift(fb);
    
    if false
        % use in special cases (very long signals)
        
        % ************************************
        % data reduction:
        % If a filter is used, details smaller than windowBW are
        % suppressed already by the filter, and using more samples
        % gives no additional information. 
        % ************************************
        if numel(fr) > p.NMax
            switch(p.filter)
              case 'gaussian'
                % 0.2 offset from the peak causes about 0.12 dB error
                incr=floor(windowBW/binspacing_Hz/5);
              case 'brickwall'
                % there is no error at all for a peak
                incr=floor(windowBW/binspacing_Hz/3);
              case 'none'
                % there is no filter, we cannot discard data at this stage
                incr=-1;
            end
            
            if incr > 1 
                fr=fr(1:incr:end);
                fb=fb(1:incr:end);
            end  
        end
    end
    
    % ************************************
    % data reduction: 
    % discard beyond fMin / fMax, if given
    % ************************************
    indexMin = 1;
    indexMax = numel(fb);  
    flag=0;
    if isfield(p, 'fMin_Hz')
        indexMin=min(find(fb >= p.fMin_Hz));
        flag=1;
    end
    if isfield(p, 'fMax_Hz')
        indexMax=max(find(fb <= p.fMax_Hz));
        flag=1;
    end
    if flag
        fb=fb(indexMin:indexMax);
        fr=fr(indexMin:indexMax);
    end
    
    if p.NMax > 0
        if p.logscale
            % ************************************
            % Sample number reduction for logarithmic
            % frequency scale
            % ************************************     
            assert(isfield(p, 'fMin_Hz'), 'need fMin_Hz in logscale mode');
            assert(p.fMin_Hz > 0, 'need fMin_Hz > 0 in logscale mode');
            if ~isfield(p, 'fMax_Hz')
                p.fMax_Hz = p.rate_Hz / 2;
            end
            
            % averaged output arrays
            fbOut=zeros(1, p.NMax-1);
            frOut=zeros(1, p.NMax-1);

            % candidate output frequencies (it's not certain yet
            % that there is actually data)
            s=logspace(log10(p.fMin_Hz), log10(p.fMax_Hz), p.NMax);

            f1=s(1);
            nextStartIndex=max(find(fb < f1));
            if isempty(nextStartIndex)
                nextStartIndex=1;
            end
            
            % iterate through all frequency regions
            % collect data
            % average
            for index=2:size(s, 2)
                f2=s(index);
                endIndex=max(find(fb < f2));

                % number of data points in bin
                n=endIndex-nextStartIndex+1;             
                if n > 0
                    % average
                    ix=nextStartIndex:endIndex;
                    fbOut(index-1)=sum(fb(ix))/n;
                    frOut(index-1)=sum(fr(ix))/n;
                    nextStartIndex=endIndex+1;
                else
                    % mark point as invalid (no data)
                    fbOut(index-1)=nan;
                end
            end
            % remove bins where no data point contributed
            ix=find(~isnan(fbOut));
            fbOut=fbOut(ix);
            frOut=frOut(ix);
            fb=fbOut;
            fr=frOut;
        else
            % ************************************
            % Sample number reduction for linear
            % frequency scale
            % ************************************     
            len=size(fb, 2);
            decim=ceil(len/p.NMax); % one sample overlength => decim=2                              
            if decim > 1
                % truncate to integer multiple
                len=floor(len / decim)*decim;
                fr=fr(1:len);
                fb=fb(1:len);

                fr=reshape(fr, [decim, len/decim]);
                fb=reshape(fb, [decim, len/decim]);
                if singleBinMode
                    % apply peak hold over each segment (column)
                    fr=max(fr);
                else
                    % apply averaging over each segment (column)
                    fr = sum(fr) / decim;
                end
                fb=sum(fb)/decim; % for each column the average
            end % if sample reduction necessary
        end % if linear scale
    end % if sample number reduction

    % ************************************
    % convert magnitude to dB.
    % truncate very small values 
    % using an artificial noise floor
    % ************************************     
    fr=(10*log10(fr+noisefloor));

    if singleBinMode
        % ************************************
        % The power reading shows the content of each 
        % FFT bin. This is accurate for single-frequency
        % tones that fall exactly on the frequency grid
        % (an integer number of cycles over the signal length)
        % ************************************     
    else
        % ************************************
        % convert sensed power density from FFT bin spacing
        % to resolution bandwidth
        % ************************************     
        fr=fr+10*log10(p.RBW_power_Hz/binspacing_Hz);
    end

    % ************************************
    % Post-processing: 
    % Translate frequency axis to RF
    % ************************************
    fb = fb + p.originMapsTo_Hz;
    
    % ************************************
    % convert to requested units
    % ************************************
    switch(p.freqUnit)
      case 'Hz'
      case 'kHz'
        fb = fb / 1e3;
      case 'MHz'
        fb = fb / 1e6;
      case 'GHz'
        fb = fb / 1e9;
      otherwise
        error('unsupported frequency unit');
    end
    
    % ************************************
    % Plot (if requested)
    % ************************************    
    if p.fig > 0
        % *************************************************************
        % title
        % *************************************************************
        if isfield(p, 'title')
            t=['"', p.title, '";'];
        else
            t='';
        end
        switch(p.filter)
          case 'gaussian'
            t=[t, ' filter: Gaussian '];
          case 'brickwall'
            t=[t, ' filter: ideal bandpass '];
          case 'none'
            t=[t, ' filter: none '];
          otherwise
            assert(0)
        end
        
        if ~strcmp(p.filter, 'none')
            t=[t, '(', mat2str(p.RBW_window_Hz), ' Hz)'];       
        end
        
        if isfield(p, 'pNom_dBm')
            % *************************************************************
            % show dB relative to a given absolute power in dBm 
            % *************************************************************
            fr=fr-p.pNom_dBm;
            yUnit='dB';
        else
            yUnit='dBm';
        end
        
        % *************************************************************
        % plot
        % *************************************************************
        figure(p.fig);    
        if p.logscale
            handle = semilogx(fb, fr, p.plotStyle);
        else
            handle = plot(fb, fr, p.plotStyle);
        end
        hold on; % after plot, otherwise prevents logscale
        
        if isfield(p, 'lineWidth')
            set(handle, 'lineWidth', p.lineWidth);
        end
        
        % *************************************************************
        % axis labels
        % *************************************************************
        xlabel(p.freqUnit);
        if singleBinMode
            ylabel([yUnit, ' (continuous wave)']);    
        else
            ylabel([yUnit, ' in ', mat2str(p.RBW_power_Hz), ' Hz integration BW']);
        end
        title(t);
        
        % *************************************************************
        % adapt y range, if requested
        % *************************************************************
        y=ylim();
        y1=y(1); y2=y(2);

        rescale=false;
        if isfield(p, 'yMin')
            y(1)=p.yMin;
            rescale=true;
        end
        if isfield(p, 'yMax')
            y(2)=p.yMax;
            rescale=true;
        end
        if rescale
            ylim(y);
        end
    end
end

function [winLeft, winRight] = spectrumAnalyzerGetWindow(filter, singleBinMode, RBW_window_Hz, binspacing_Hz)
    
    switch(filter)
      case 'gaussian'
        
        % construct Gaussian filter
        % -60 / -3 dB shape factor 4.472
        nRBW=6;
        nOneSide=ceil(RBW_window_Hz/binspacing_Hz*nRBW); 
        
        filterBase=linspace(0, nRBW, nOneSide);
        winLeft=exp(-(filterBase*0.831) .^2);
        winRight=fliplr(winLeft(2:end)); % omit 0 Hz bin
        
      case 'brickwall'
        nRBW=1;
        n=ceil(RBW_window_Hz/binspacing_Hz*nRBW); 
        n1 = floor(n/2);
        n2 = n - n1;
        
        winLeft=ones(1, n1);
        winRight=ones(1, n2); 
        
      case 'none'
        winLeft=[];
        winRight=[];
      
      otherwise
        error('unknown RBW filter type');
    end
    
    % the window is not supposed to shift the spectrum.
    % Therefore, the first bin is always in winLeft (0 Hz):
    if size(winLeft, 2)==1 && isempty(winRight)
        % there is no use to convolve with one-sample window
        % it's always unity
        winLeft=[];
        winRight=[];
        tmpwin=[];
    end
    
    if ~isempty(winLeft)
        % (note: it is not possible that winRight is empty, while winLeft is not)
        if singleBinMode
            % normalize to unity at 0 Hz
            s=winLeft(1);
        else
            % normalize to unity area under the filter
            s=sum(winLeft)+sum(winRight);
        end
        winLeft=winLeft / s;
        winRight=winRight / s;
    end
end

% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
%   varargin = {'hello'}
    r = flattenCellArray(varargin, struct()); 
end

function r = flattenCellArray(arr, r)
    ix=1;
    ixMax = numel(arr);
    while ix <= ixMax
        e = arr{ix};
        
        if iscell(e)
            % cell array at 'key' position gets recursively flattened 
            % becomes struct
            r = flattenCellArray(e, r);
        elseif ischar(e)
            % string => key. 
            % The following entry is a value
            ix = ix + 1;
            v = arr{ix};
            % store key-value pair
            r.(e) = v;
        elseif isstruct(e)
            names = fieldnames(e);
            for ix2 = 1:numel(names) 
                k = names{ix2};
                r.(k) = e.(k);
            end
        else 
            e
            assert(false)
        end
        ix=ix+1;
    end % while
end