## 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

%add up subfilters, advancing f3 one stage
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
%                   'ill'        ILL (Adobe Illustrator)
%                   '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;

m = Theta_LS(1);
b = Theta_LS(2);

## Measuring Peak-to-Peak

September 30, 2011 Coded in Matlab
function y = peak2peak(signal)
% Return the peak-to-peak amplitude of the supplied signal.  This is the
% same as max(signal) minus min(signal).
%
% Usage: y = PEAK2PEAK(signal);
%
%        SIGNAL is your one-dimensional input array
%
% Author: sparafucile17

% Input must have some length
if(length(signal) == 1)
error('ERROR: input signal must have more than one element');
end

% This function only supports one-dimensional arrays
if((size(signal, 2) ~= 1) && (size(signal, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end

% Find the peak and return it
min_sig = min(signal);
max_sig = max(signal);

y = max_sig - min_sig;

## Percent Error

September 30, 2011 Coded in Matlab
function [per_error] = percent_error(measured, actual)
% Creates the percent error between measured value and actual value
%
% Usage: percent_error(MEASURED, ACTUAL);
%
%        Measured is the your result
%        Actual is the value that your result should be
%
% Author: sparafucile17

per_error = abs(( (measured - actual) ./ actual ) * 100);

## Finding local minima

September 26, 20111 comment Coded in Matlab
function [time,amp] = minima(signal, len)
%
% Finds the local minimum values of the given signal.
%
% Usage:  [IND,AMP] = MINIMA(X, N);
%
%         N is the number of points that need to be evaluated
%           (Normally equal to LENGTH(X))
%         X is the data
%         IND contains the indices of X that contain the minima data
%         AMP contains the minima amplitudes
%
% Author: sparafucile17 10/02/2003

%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1;  %allow a maxima point at the second sample

%prevent length error
if(len > length(signal))
len = length(signal)
end

%Loop through data and find local minimums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0))  %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end

%After loop assign data to output variables
time = local_time;
amp = local_amp;

## Finding local maxima

September 26, 20112 comments Coded in Matlab
function [time,amp] = maxima(signal, thresh, len)
% Finds the local maximum values of the given signal.
%
% Usage:  [IND,AMP] = MAXIMA(X, THRESH, N);
%
%         X is the data
%         THRESH is the threshold that signal must be greater
%            than in order to be counted. (Default = 0.0)
%         N is the number of points that need to be evaluated
%           (Defaults = LENGTH(X))
%         IND contains the indices of X that contain the maxima data
%         AMP contains the maxima amplitudes
%
% Author: sparafucile17 10/02/2003

if(exist('len') == 0)
len = length(signal);
end
if(exist('thresh') == 0)
thresh = 0.0;
end

%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1;  %allow a maxima point at the second sample

%prevent length error
if(len > length(signal))
len = length(signal)
end

%Loop through data and find local maximums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0) & (cur > thresh))  %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end

%After loop assign data to output variables
time = local_time;
amp = local_amp;