## Computing Optimum Two-Stage Decimation Factors

November 21, 2011 Coded in Matlab
``````function [D1,D2] = Dec_OptimTwoStage_Factors(D_Total,Passband_width,Fstop)
%
%    When breaking a single large decimation factor (D_Total) into
%    two stages of decimation (to minimize the orders of the
%    necessary lowpass digital filtering), an optimum first
%    stage of decimation factor (D1) can be found.  The second
%    stage's decimation factor is then D_Total/D1.
%
%    Inputs:
%         D_Total = original single-stage decimation factor.
%         Passband_width = desired passband width of original
%                          single-stage lowpass filter (Hz).
%         Fstop = desired beginning of stopband freq of original
%                          single-stage lowpass filter (Hz).
%
%    Outputs:
%         D1 = optimum first-stage decimation factor.
%         D2 = optimum second-stage decimation factor.
%
%    Example: Assume you want to decimate a sequence by D_Total=100,
%             your original single-stage lowpass filter's passband
%             width and stopband frequency are 1800 and 2200 Hz
%             respectively.  We use:
%
%              [D1,D2] = Dec_OptimTwoStage_Factors(100, 1800, 2200)
%
%             giving us the desired D1 = 25, and D2 = 4.
%             (That is, D1*D2 = 25*4 = 100 = D_Total.)
%
%    Author: Richard Lyons, November, 2011

% Start of processing
DeltaF = (Fstop -Passband_width)/Fstop;

Radical = (D_Total*DeltaF)/(2 -DeltaF); % Terms under sqrt sign.
Numer = 2*D_Total*(1 -sqrt(Radical));  % Numerator.
Denom = 2 -DeltaF*(D_Total + 1);  % Denominator.
D1_estimated = Numer/Denom;  %Optimum D1 factor, but not
%                  an integer.

% Find all factors of 'D_Total'
Factors_of_D_Total = Find_Factors_of_D_Total(D_Total);

% Find the one factor of 'D_Total' that's closest to 'D1_estimated'
Temp = abs(Factors_of_D_Total -D1_estimated);
Index_of_Min = find(Temp == min(Temp)); % Index of optim. D1
D1 = Factors_of_D_Total(Index_of_Min); % Optimum first
%    decimation factor
D2 = D_Total/D1;  % Second decimation factor.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allfactors] = Find_Factors_of_D_Total(X)
% Finds all integer factors of intger X.
% Filename was originally 'listfactors.m', written
[b,m,n] = unique(factor(X));
%b is all prime factors
occurences = [m(1) diff(m)];
current_factors = [b(1).^(0:occurences(1))]';
for index = 2:length(b)
currentprime = b(index).^(0:occurences(index));
current_factors = current_factors*currentprime;
current_factors = current_factors(:);
end
allfactors = sort(current_factors);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%``````

## FFT for arbitrary time and frequency grids (Chirp-z Transform)

November 17, 2011 Coded in Matlab
``````function Demo

while 1

%Frequency grid with arbitrary start value, spacing, and number of elements. ...
%Mathematically, the grid is
%
%   fo+Df*m, m=0, 1,...,M-1.
%

fo = 10*rand(1)-0.5; %Start value.
Df = rand(1)+0.01;   %Frequency spacing.
M = round(rand(1)*1000); %Number of frequencies.

%Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
%the grid is
%
%   to+Dt*n, n=0, 1,...,N-1.
%

to = 10*rand(1)-0.5; %Start value; (instant of first sample).
Dt = rand(1)+0.01;   %Time spacing.
N = round(rand(1)*1000); %Number of samples.

%Random vector of samples.
x = randn(1,N)+1i*randn(1,N);

%x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N.

%We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.

%Compute DFT using the direct and chirp-z methods.

tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;

disp(['Timing direct method: ', num2str(tmD), ' sec.'])
disp(['Timing chirp z: ',num2str(tmF),' sec.'])
disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
disp(' ')
input('(Press a key for another trial or Ctrl-C) ')
disp(' ')
end

function X = chirpzDFT(x,to,Dt,fo,Df,M)

%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
%               N
%       X(k) = sum  x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
%              n=1
%
%Input parameters:
%
%  x: Signal samples.
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
%  M: Desired number of output samples.
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%

x = x(:).';

N = numel(x);
P = 2^ceil(log2(M+N-1));

%The next three lines do not depend on the vector x, and so they can be pre-computed if
%the time and frequency grids do not change, (i.e, for computing the transform of
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and
%three if the grids change.

a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse.

%Weigh x using a and perform FFT convolution with phi.
X = ifft( fft(x.*a,P) .* phi );

%Truncate the convolution tails and weigh using b.
X = X(N:M+N-1) .* b;

function X = DirectGridDFT(x,to,Dt,fo,Df,M)

%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.

x = x(:).';

N = length(x);

X = zeros(1,M);

for k = 1:M
for n = 1:N
X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
end
end``````

## NCO bit true model in a single equation

October 23, 2011 Coded in Matlab
``````%NCO bit true model
clear all; close all;

%example NCO parameters
n = 20000;                     %number of test samples
A = 2^15 - 1;                  %max amplitude, amplitude resolution 16 bit signed
M = 2^32;                      %NCO accumulator phase resolution
inc = round(M*.001);           %NCO accumulator phase increment
k = 2^12;                      %lut phase resolution (lut size), may equal M

%single equation for NCO bit true model
y1 = round(A*sin(2*pi*round(k*mod((0:n-1)*inc/M,M))/k));

%%%%%%%%%%%%%%%%%%%%%%%%% equivalent code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lut = round(A*sin(2*pi*(0:k-1)/k));      %lut, one cycle sine data

ptr = 0;
for i = 1:n
ptr = mod(ptr + inc, M);      %phase accumulator(ramp)
end

%compare
plot(y1);hold;
plot(y2,'r--')
plot(y1 - y2,'g--')``````

## Complex filter with three multipliers per tap

October 22, 2011 Coded in Matlab
``````%random input
I_in = randn(1,1024);
Q_in = randn(1,1024);

%random test filter
h = randn(1,20)+j*randn(1,20);

%split up into 3 subfilters
h2 = real(h);
h3 = real(h)-imag(h);         %extra subtractor

%apply subfilter
ch1 = filter(h1,1,I_in);
ch2 = filter(h2,1,Q_in-I_in); %extra subtractor
ch3 = filter(h3,1,Q_in);

%combine outputs
yr = round((ch3-ch2)*2^15);
yi = round((ch1+ch2)*2^15);
IQ_out = complex(yr,yi);

%direct computation
ref = round(2^15*filter(h,1,complex(I_in,Q_in)));
plot(IQ_out - ref);``````

## Even/Odd FIR filter structure

October 22, 20111 comment Coded in Matlab
``````% odd/even filter structure verification
% two filter structures compared
% (1) direct use of filter function
% (2) even/odd split structure

x = randn(1,2048);
h = fir1(23,.3);
y1 = filter(h,1,x); % direct filter for reference

he = h(1:2:end); % even coeff, relative to hardware definition
ho = h(2:2:end); % odd coeffs
xe = x(1:2:end); % even samples
xo = x(2:2:end); % odd samples

f1 = filter(ho,1,xe);   %subfilter 1
f2 = filter(he,1,xo);   %subfilter 2
f3 = filter(he,1,xe);   %subfilter 3
f4 = filter(ho,1,xo);   %subfilter 4

ye = [0 f1] + [0 f2];
yo = [f3 0] + [0 f4];

%interleave even/odd into one stream
y2 = [y1 0 0];
y2(1:2:end) = ye;
y2(2:2:end) = yo;

%quantise and compare both ref and subfilters
y1 = round(y1*2^15);
y2 = round(y2*2^15);
plot(y1 - y2(2:end-1))``````

## Vector alignment for verification

October 22, 2011 Coded in Matlab
``````%"align_run" is top level for testing "align" function

function align_run
for i = 1:1
L1 = 3400;   %length of first vector
L2 = 4023;   %length of second vector

%generate two random vectors
x1 = randn(1,L1);
x2 = randn(1,L2);

%generate random index values for two identical segments
r1 = round(rand(1)*L1);
r2 = round(rand(1)*L2);
r1(r1 == 0) = 1;
r2(r2 == 0) = 1;

%length of identical segments as percent of vector length
n = round(50/100 * min([L1 L2]));
n(n >= min([L1 L2])) = min([L1 L2])-1;
while r1+n > L1, r1 = r1-1; end;
while r2+n > L2, r2 = r2-1; end;

%create identical segments at random location
x1(r1:r1+n) = x2(r2:r2+n);

%call align function
[y1,y2] = align(x1,x2);

%check alignment
subplot(2,1,1);
plot(y1);hold
plot(y2,'r--') ;
legend('y1','y2')

subplot(2,1,2);
plot(y1 - y2,'g')
legend('(y1 - y2)')

end

return

%%%%%%%%%%%%%%%%% align %%%%%%%%%%%%%%%%%%
function [y1,y2] = align(x1,x2)

%force a column orientation
x1 = x1(:);
x2 = x2(:);

L1 = length(x1);
L2 = length(x2);
L = max([L1 L2]);

%equalise vector lengths
if L1 > L2, x2 = [x2; zeros(L1-L2,1)]; end;
if L1 < L2, x1 = [x1; zeros(L2-L1,1)]; end;

%apply correlation and find index of maximum output
y = xcorr(x1,x2);
n = find(y == max(y)) ;
n = mod(n,L)+1;

%align x1 vector with x2
y1 = [x1(n:end); x1(1:n-1)] ;
y2 = x2;

return``````

## integrate RMS phase error from spectrum

October 21, 20111 comment Coded in Matlab
``````% ****************************************************************
% RMS phase error from phase noise spectrum
% Markus Nentwig, 2011
%
% - integrates RMS phase error from a phase noise power spectrum
% - generates an example signal and determines the RMS-average
%   phase error (alternative method)
% ****************************************************************
function pn_snippet()
close all;

% ****************************************************************
% Phase noise spectrum model
% --------------------------
% Notes:
% - Generally, a phase noise spectrum tends to follow
%   |H(f)| = c0 + c1 f^n1 + c2 f^n2 + c3 f^n3 + ...
%   Therefore, linear interpolation should be performed in dB
%   on a LOGARITHMIC frequency axis.
% - Single-/double sideband definition:
%   The phase noise model is defined for -inf <= f <= inf
%   in other words, it contributes the given noise density at
%   both positive AND negative frequencies.
%   Assuming a symmetrical PN spectrum, the model is evaluated
%   for |f| => no need to explicitly write out the negative frequency
%   side.
% ****************************************************************

% PN model
% first column:
%   frequency offset
% second column:
%   spectral phase noise density, dB relative to carrier in a 1 Hz
%   observation bandwidth

d = [0 -80
10e3 -80
1e6 -140
9e9 -140];
f_Hz = d(:, 1) .';
g_dBc1Hz = d(:, 2) .' -3;

% get RMS phase error by integrating the power spectrum
% (alternative 1)
e1_degRMS = integratePN_degRMS(f_Hz, g_dBc1Hz)

% get RMS phase error based on a test signal
% (alternative 2)
n = 2 ^ 20;
deltaF_Hz = 2;
e2_degRMS = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz)
end

% ****************************************************************
% Integrate the RMS phase error from the power spectrum
% ****************************************************************
function r = integratePN_degRMS(f1_Hz, g1_dBc1Hz)
1;
% integration step size
deltaF_Hz = 100;

% list of integration frequencies
f_Hz = deltaF_Hz:deltaF_Hz:5e6;

% interpolate spectrum on logarithmic frequency, dB scale
% unit is dBc in 1 Hz, relative to a unity carrier
fr_dB = interp1(log(f1_Hz+eps), g1_dBc1Hz, log(f_Hz+eps), 'linear');

% convert to power in 1 Hz, relative to a unity carrier
fr_pwr = 10 .^ (fr_dB/10);

% scale to integration step size
% unit: power in deltaF_Hz, relative to unity carrier
fr_pwr = fr_pwr * deltaF_Hz;

% evaluation frequencies are positive only
% phase noise is two-sided
% (note the IEEE definition: one-half the double sideband PSD)
fr_pwr = fr_pwr * 2;

% sum up relative power over all frequencies
pow_relToCarrier = sum(fr_pwr);

% convert the phase noise power to an RMS magnitude, relative to the carrier
pnMagnitude = sqrt(pow_relToCarrier);

% convert from radians to degrees
r = pnMagnitude * 180 / pi;
end

% ****************************************************************
% Generate a PN signal with the given power spectrum and determine
% the RMS phase error
% ****************************************************************
function r = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz);
A = 1; % unity amplitude, arbitrary

% indices of positive-frequency FFT bins.
% Does not include DC
ixPos = 2:(n/2);

% indices of negative-frequency FFT bins.
% Does not include the bin on the Nyquist limit
assert(mod(n, 2) == 0, 'n must be even');
ixNeg = n - ixPos + 2;

% Generate signal in the frequency domain
sig = zeros(1, n);
% positive frequencies
sig(ixPos) = randn(size(ixPos)) + 1i * randn(size(ixPos));
% for purely real-valued signal: conj()
% for purely imaginary-valued signal such as phase noise: -conj()
% Note:
%   Small-signals are assumed. For higher "modulation indices",
%   phase noise would become a circle in the complex plane
sig(ixNeg) = -conj(sig(ixPos));

% normalize energy to unity amplitude A per bin
sig = sig * A / RMS(sig);

% normalize energy to 0 dBc in 1 Hz
sig = sig * sqrt(deltaF_Hz);

% frequency vector corresponding to the FFT bins (0, positive, negative)
fb_Hz = FFT_freqbase(n, deltaF_Hz);

% interpolate phase noise spectrum on frequency vector
% interpolate dB on logarithmic frequency axis
H_dB = interp1(log(f_Hz+eps), g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');

% dB => magnitude
H = 10 .^ (H_dB / 20);

% plot on linear f axis
figure(); hold on;
plot(fftshift(fb_Hz), fftshift(mag2dB(H)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (linear frequency axis)');

% plot on logarithmic f axis
figure(); hold on;
ix = find(fb_Hz > 0);
semilogx(fb_Hz(ix), mag2dB(H(ix)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (logarithmic frequency axis)');

% apply frequency response to signal
sig = sig .* H;

sig (1) = A;

% convert to time domain
td = ifft(sig);

% determine phase
% for an ideal carrier, it would be zero
% any deviation from zero is phase error
ph_deg = angle(td) * 180 / pi;

figure();
ix = 1:1000;
plot(ix / n / deltaF_Hz, ph_deg(ix));
xlabel('time / s');
ylabel('phase error / degrees');
title('phase of example signal (first 1000 samples)');

figure();
hist(ph_deg, 300);
title('histogram of example signal phase');
xlabel('phase error / degrees');

% RMS average
% note: exp(1i*x) ~ x
%       as with magnitude, power/energy is proportional to phase error squared
r = RMS(ph_deg);

% knowing that the signal does not have an amplitude component, we can alternatively
% determine the RMS phase error from the power of the phase noise component

% exclude carrier:
sig(1) = 0;

% 3rd alternative to calculate RMS phase error
rAlt_deg = sqrt(sig * sig') * 180 / pi
end

% gets 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

% Root-mean-square average
function r = RMS(vec)
r = sqrt(vec * vec' / numel(vec));
end

% convert a magnitude to dB
function r = mag2dB(vec)
r = 20*log10(abs(vec) + eps);
end``````

## Saving a Matlab figure to file using command-line

September 30, 2011 Coded in Matlab
``````function save_fig(h, name, format)
% Matlab has a wierd way of saving figures and keeping the plot looking
% exactly the way it does on your screen.  Thi function simplifies the
% entire process in a way that is barely documented.
%
% Usage:      SAVE_FIG(h, name, format);
%
%             H is the handle to the figure.  This can be obtain in the
%               following manner:  H = figure(1);
%             NAME is the filename without extension
%             FORMAT is the graphic format.  Options are:
%
%                   'bmpmono'    BMP monochrome BMP
%                   'bmp16m'     BMP 24-bit BMP
%                   'bmp256'     BMP 8-bit (256-color)
%                   'bmp'        BMP 24-bit
%                   'meta'       EMF
%                   'eps'        EPS black and white
%                   'epsc'       EPS color
%                   'eps2'       EPS Level 2 black and white
%                   'epsc2'      EPS Level 2 color
%                   'fig'        FIG Matlab figure file
%                   'hdf'        HDF 24-bit
%                   'jpeg'       JPEG 24-bit
%                   'pbm'        PBM (plain format) 1-bit
%                   'pbmraw'     PBM (raw format) 1-bit
%                   'pcxmono'    PCX 1-bit
%                   'pcx24b'     PCX 24-bit color PCX file format
%                   'pcx256'     PCX 8-bit newer color (256-color)
%                   'pcx16'      PCX 4-bit older color (16-color)
%                   'pdf'        PDF Color PDF file format
%                   'pgm'        PGM Portable Graymap (plain format)
%                   'pgmraw'     PGM Portable Graymap (raw format)
%                   'png'        PNG 24-bit
%                   'ppm'        PPM Portable Pixmap (plain format)
%                   'ppmraw'     PPM Portable Pixmap (raw format)
%                   'svg'        SVG Scalable Vector Graphics
%                   'tiff'       TIFF 24-bit
%
% Author: sparafucile17

if(strcmp(format,'fig'))
saveas(h, name, 'fig');
else
options.Format = format;
hgexport(h, name, options);
end``````

## Putting a Matlab figure to a specific screen location

September 30, 2011 Coded in Matlab
``````function set_fig_position(h, top, left, height, width)
% Matlab has a wierd way of positioning figures so this function
% simplifies the poisitioning scheme in a more conventional way.
%
% Usage:      SET_FIG_POISITION(h, top, left, height, width);
%
%             H is the handle to the figure.  This can be obtain in the
%               following manner:  H = figure(1);
%             TOP is the "y" screen coordinate for the top of the figure
%             LEFT is the "x" screen coordinate for the left side of the figure
%             HEIGHT is how tall you want the figure
%             WIDTH is how wide you want the figure
%
% Author: sparafucile17

% PC's active screen size
screen_size = get(0,'ScreenSize');
pc_width  = screen_size(3);
pc_height = screen_size(4);

%Matlab also does not consider the height of the figure's toolbar...
%Or the width of the border... they only care about the contents!
toolbar_height = 77;
window_border  = 5;

% The Format of Matlab is this:
%[left, bottom, width, height]
m_left   = left + window_border;
m_bottom = pc_height - height - top - toolbar_height - 1;
m_height = height;
m_width  = width - 1;

%Set the correct position of the figure
set(h, 'Position', [m_left, m_bottom, m_width, m_height]);

%If you want it to print to file correctly, this must also be set
% and you must use the "-r72" scaling to get the proper format
set(h, 'PaperUnits', 'points');
set(h, 'PaperSize', [width, height]);
set(h, 'PaperPosition', [0, 0, width, height]); %[ left, bottom, width, height]``````

## Least Squares fit to a linear equation

September 30, 2011 Coded in Matlab
``````function [m, b] = ls_linear(x_observed, y_observed)
%
% Perform Least Squares curve-fitting techniques to solve the coefficients
% of the linear equation: y = m*x + b.  Input and Output equation
% observations must be fed into the algorithm and a best fit equation will
% be calculated.
%
% Usage:     [M,B] = ls_linear(observed_input, observed_output);
%
%            observed_input  is the x-vector observed data
%            observed_output is the y-vector observed data
%
% Author:    sparafucile17 06/04/03
%

if(length(x_observed) ~= length(y_observed))
error('ERROR: Both X and Y vectors must be the same size!');
end

%Calculate vector length once
vector_length = length(x_observed);

%
% Theta = [   x1    1   ]
%         [   x2    1   ]
%         [   x3    1   ]
%         [   x4    1   ]
%         [   ...  ...  ]
%
Theta = ones(vector_length, 2);
Theta(1:end, 1) = x_observed;

%
% Theta_0 = [  y1   ]
%           [  y2   ]
%           [  y3   ]
%           [  y4   ]
%           [  ...  ]
%
Theta_0 = y_observed;

%
% Theta_LS = [  M  ]
%            [  B  ]
%
Theta_LS = ((Theta' * Theta)^-1 ) * Theta' * Theta_0;