DSPRelated.com
   Code Snippets

Even/Odd FIR filter structure

kaz - ● October 22, 2011●1 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

kaz - ● 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

Markus Nentwig ● October 21, 2011●1 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;
    
    % add carrier
    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

Jeff T ● 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

Jeff T ● 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

Jeff T ● 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;

%Return the answer
m = Theta_LS(1);
b = Theta_LS(2);

Measuring Peak-to-Peak

Jeff T ● 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);

% Answer
y = max_sig - min_sig;

Percent Error

Jeff T ● 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

Jeff T ● September 26, 2011●1 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

Jeff T ● September 26, 2011●2 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;