DSP Code Snippets

● November 2, 2010 ● Coded in Matlab

%==============================================================================
%
%                   Gray Code Generator for Verilog
%
%   This script generates a Verilog header (.vh) file which contains a set of
%   of `define directives to establish Gray-coded states for a state machine.
%   The output of this script is intended to be used in conjunction with a
%   Verilog design which uses Gray coding in some application.  For more
%   information on Gray coding, please see 
%   http://mathworld.wolfram.com/GrayCode.html.
%
%   Usage: 
%
%       The bit width of the Gray coded states is set as a script variable.
%       2^bit_width states are generated.  The output filename is also 
%       set as a script variable.
%
%   Author: Tom Chatt
%
%------------------------------------------------------------------------------
% Revision History:
%
%   Original Version        Tom Chatt           11/2/2010
%
%==============================================================================

clear all; clc;

%------------------------------------------------------------------------------

bit_width = 8;  % Bit width of the generated code

% Print the generated code entries to a Verilog defines file?
print_to_v_def_file = true;

output_file = ['gray_code_',num2str(bit_width),'.vh']; % Output file name

% Gray code state name prefix.  State names will be of the form 
% GC<bit_width>B_<state number>.
prefix = ['GC',num2str(bit_width),'B_'];

%--------------Do not edit code below this point-------------------------------

% Generate the code, starting with 1 bit code and looping upward
code = [0; 1];
if bit_width ~= 1
    for i = 2 : bit_width
        code = [zeros(size(code,1), 1), code; ones(size(code,1), 1), flipud(code)]; %#ok<AGROW>
    end
end

% Print out the generated code
fid = fopen(output_file, 'w');

for i = 1 : size(code,1)
   
    macro_name = [prefix, num2str(i)];
    macro_val = [num2str(bit_width),'''', dec2bin(binvec2dec(code(i,:)),bit_width)];
    fprintf(fid, ['ifndef ',macro_name,'\n']);
    fprintf(fid, ['    ','`define ', macro_name, ' ', macro_val, '\n']);
    fprintf(fid, 'endif\n');
    
end

fclose('all');

DWT frequency division graphing

David Valencia ● October 29, 2010●4 comments ● Coded in Matlab

% ----------------------------------------------------------
% Title: Discrete Wavelet Transform
% Author: David Valencia
% UPIITA IPN 2010
%
% For DSPRelated.com
% in: http://www.dsprelated.com/showcode/13.php
%
% Computes the Discrete Wavelet Transform
% equivalent filters
%
% Dependencies: 
% formfilter.m - http://www.dsprelated.com/showcode/12.php
% upsample2.m - http://www.dsprelated.com/showcode/10.php
%
% Revisions:
%       v1.0a Commented and translated in English
%
% More information at:
% http://www.dsprelated.com/showarticle/115.php
% -----------------------------------------------------------

close all; clear all; clc;

disp('== GENERALIZED DWT FILTERS ==')

%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db6';

% Obtain base filters
if(typeofbasis == 'b')
    [Rf,Df] = biorwavf(typbior);
    [h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
    [h0,h1,g0,g1] = wfilters(typor);
end;

%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.

L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches

% L = Basic filters length (h0 ó h1)
% N = Input vector length
% n_stages = stage quantity, it generates (2^n_stages) branches

figure;
[H1,Wg]=freqz(h1,1,512);
plot(Wg,abs(H1),'r','linewidth',2); hold on;

for i = 0:(n_branches/2)-1
    if n_stages < 3
        hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
    else
        hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
    end
    Lhx = length(hx);
    
    % Graphing code
    [H0,Wg]=freqz(hx,1,512);
    plot(Wg,abs(H0),'b','linewidth',2); hold on;
    pause; % Used to see how each filter appears
end

checking resampling in time domain

kaz - ● January 21, 2012 ● Coded in Matlab

clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,512));
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor

%%%%%%%%%%%%%%%%%%%%% up/dn model %%%%%%%%%%%%%%%%%%%%%%
%upsample input
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f = filter(up*h,1,x_up);

%downsample signal by decimation
x_dn = x_f(1:dn:end);

delay = 30/2+1;  
figure;
subplot(2,1,1);hold;
plot(x_up,'o--');
plot(x_f(delay:end),'r.-');
legend('input with zeros','upsampled stage');

subplot(2,1,2);hold;
plot(x_up(1:dn:end),'o--');
plot(x_dn(ceil(delay/dn):end),'r.-');
legend('input with zeros','final signal');

Baseband-equivalent phase noise model

Markus Nentwig ● December 18, 2011●4 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

Peak-to-average ratio analysis

Markus Nentwig ● December 10, 2011 ● Coded in Matlab

% *************************************************************
% Peak-to-average analyzis of complex baseband-equivalent signal
% Markus Nentwig, 10/12/2011
% Determines the magnitude relative to the RMS average
% which is exceeded with a given (small) probability
% A typical probability value is 99.9 %, which is very loosely related
% to an uncoded bit error rate target of 0.1 %.
% *************************************************************
function sn_headroom()
    % number of samples for example signals
    n = 1e6;

    % *************************************************************
    % example: 99.9% PAR for white Gaussian noise
    % *************************************************************
    signal = 12345 * (randn(1, n) + 1i * randn(1, n));
    [PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, false);
    printf('white Gaussian noise: %1.1f dB PAR at radio frequency\n', PAR_dB);
    [PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, true);
    printf('white Gaussian noise: %1.1f dB PAR at baseband\n', PAR_dB);

    % *************************************************************
    % example: PAR of continuous-wave signal
    % *************************************************************
    % a quarter cycle gives the best coverage of amplitude values
    % the statistics of the remaining three quarters are identical (symmetry)
    signal = 12345 * exp(1i*2*pi*(0:n-1) / n / 4);

    % Note: peaks can be below the average, depending on the given probability
    % a possible peak-to-average ratio < 1 (< 0 dB) is a consequence of the
    % peak definition and not unphysical.
    % For a continuous-wave signal, the average value equals the peak value, 
    % and PAR results slightly below 0 dB are to be expected.
    [PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, false);
    printf('continuous-wave: %1.1f dB PAR at radio frequency\n', PAR_dB);
    [PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, true);
    printf('continuous-wave: %1.1f dB PAR at baseband\n', PAR_dB);

    % *************************************************************
    % self test:
    % For a real-valued Gaussian random variable, the probability
    % to not exceed n sigma is
    % n = 1: 68.27 percent 
    % n = 2: 95.5 percent
    % n = 3: 99.73 percent
    % see http://en.wikipedia.org/wiki/Normal_distribution
    % section "confidence intervals"
    % *************************************************************
    signal = 12345 * randn(1, n);
    [PAR, PAR_dB] = peakToAverageRatio(signal, 68.2689/100, false);
    printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(1), PAR_dB);
    [PAR, PAR_dB] = peakToAverageRatio(signal, 95.44997/100, false);
    printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(2), PAR_dB);
    [PAR, PAR_dB] = peakToAverageRatio(signal, 99.7300/100, false);
    printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(3), PAR_dB);
end

function [PAR, PAR_dB] = peakToAverageRatio(signal, peakProbability, independentIQ)
    1;
    % force row vector
    assert(min(size(signal)) == 1, 'matrix not allowed');
    signal = signal(:) .';

    assert(peakProbability > 0);
    assert(peakProbability <= 1);
    
    % *************************************************************
    % determine RMS average and normalize signal to unity
    % *************************************************************
    nSamples = numel(signal);
    sAverage = sqrt(signal * signal' / nSamples);
    signal = signal / sAverage;

    % *************************************************************
    % "Peaks" in a complex-valued signal can be defined in two
    % different ways:
    % *************************************************************
    if ~independentIQ
        % *************************************************************
        % Radio-frequency definition: 
        % ---------------------------
        % The baseband-equivalent signal represents the modulation on 
        % a radio frequency carrier
        % The instantaneous magnitude of the modulated radio frequency 
        % wave causes clipping, for example in a radio frequency 
        % power amplifier. 
        % The -combined- magnitude of in-phase and quadrature signal 
        % (that is, real- and imaginary part) is relevant.
        % This is the default definition, when working with radio
        % frequency (or intermediate frequency) signals, as long as a 
        % single, real-valued waveform is being processed.
        % *************************************************************
        sMag = abs(signal); 
        t = 'mag(s) := |s|; cdf(mag(s))';
    else
        % *************************************************************
        % Baseband definition
        % -------------------
        % The baseband-equivalent signal is explicitly separated into
        % an in-phase and a quadrature branch that are processed 
        % independently.
        % The -individual- magnitudes of in-phase and quadrature branch
        % cause clipping. 
        % For example, the definition applies when a complex-valued
        % baseband signal is processed in a digital filter with real-
        % valued coefficients, which is implemented as two independent,
        % real-valued filters on real part and imaginary part.
        % *************************************************************
        % for each sample, use the maximum of in-phase / quadrature.
        % If both clip in the same sample, it's counted only once.
        sMag = max(abs(real(signal)), abs(imag(signal)));
        t = 'mag(s) := max(|real(s)|, |imag(s)|); cdf(mag(s))';
    end

    % *************************************************************
    % determine number of samples with magnitude below the "peak" 
    % level, based on the given peak probability
    % for example: probability = 0.5 => nBelowPeakLevel = nSamples/2
    % typically, this will be a floating point number
    % *************************************************************
    nBelowPeakLevel = peakProbability * nSamples;
    
    % *************************************************************
    % interpolate the magnitude between the two closest samples
    % *************************************************************
    sMagSorted = sort(sMag);
    x = [0 1:nSamples nSamples+1];
    y = [0 sMagSorted sMagSorted(end)];
    magAtPeakLevel = interp1(x, y, nBelowPeakLevel, 'linear');
        
    % *************************************************************
    % Peak-to-average ratio
    % signal is normalized, average is 1
    % the ratio relates to the sample magnitude
    % "power" is square of magnitude => use 20*log10 for dB conversion.
    % *************************************************************
    PAR = magAtPeakLevel;
    PAR_dB = 20*log10(PAR + eps);

    % *************************************************************
    % illustrate the CDF and the result
    % *************************************************************
    if true
        figure();
        plot(y, x / max(x)); 
        title(t);
        xlabel('normalized magnitude m');
        ylabel('prob(mag(s) <  m)');
        line([1, 1] * magAtPeakLevel, [0, 1]);
        line([0, max(y)], [1, 1] * peakProbability);
    end
end

Convolutional coding using transform domain approach

Senthilkumar ● March 28, 2012 ● Coded in Scilab

function [x1D,x2D]= ConvolutionCode_TransDomain()         
//Convolutional code - Transform domain approach
//g1D =  generator polynomial 1
//g2D =  generator polynomial 2
//x1D = top output sequence polynomial
//x2D = bottom output sequence polynomial 
D = poly(0,'D');
g1D = input('Enter the generator polynomial 1=') //generator polynomial 1
g2D = input('Enter the generator polynomial 2=') //generator polynomial 2
mD = input('Enter the message sequence')//message sequence polynomial representation
x1D = g1D*mD; //top output polynomial
x2D = g2D*mD; //bottom output polynomial
x1 = coeff(x1D);
x2 = coeff(x2D);
disp(modulo(x1,2),'top output sequence')
disp(modulo(x2,2),'bottom output sequence')
endfunction
//Result
//Enter the generator polynomial 1 =  1+D+D^2
//Enter the generator polynomial 2 =  1+D^2;
//Enter the message sequence = 1+0+0+D^3+D^4;
//top output sequence   
// 
//    1.    1.    1.    1.    0.    0.    1.   
//bottom output sequence   
// 
//    1.    0.    1.    1.    1.    1.    1.  
//x2D  =
// 
//         2   3   4   5   6  
//    1 + D + D + D + D + D   
//x1D  =
// 
//             2   3    4    5   6  
//    1 + D + D + D + 2D + 2D + D

Auto Wah audio effect (LFO control)

Gabriel Rivas ● September 24, 2011 ● Coded in C

/*************** AutoWah.c ********************************/

#include "bp_iir.h"
#include "AutoWah.h"

static short center_freq;
static short samp_freq;
static short counter;
static short counter_limit;
static short control;
static short max_freq;
static short min_freq;
static short f_step;
static struct bp_filter H;

/*
This is the auto wah effect initialization function. 
This initializes the band pass filter and the effect control variables
*/
void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step) {
	double C;

	/*Process variables*/
	center_freq = 0;
	samp_freq = sampling;
	counter = effect_rate;
	control = 0;

	/*User Parametters*/
	counter_limit = effect_rate;
	
	/*Convert frequencies to index ranges*/
	min_freq = 0;
	max_freq = (maxf - minf)/freq_step;

        bp_iir_init(sampling,gainfactor,Q,freq_step,minf);
	f_step = freq_step; 
}

/*
This function generates the current output value
Note that if the input and output signal are integer
unsigned types, we need to add a half scale offset
*/
double AutoWah_process(int xin) {
	double yout;

	yout = bp_iir_filter(xin,&H);
	#ifdef INPUT_UNSIGNED
		yout += 32767;
	#endif
	
	return yout;
}

/*
This function will emulate a LFO that will vary according
to the effect_rate parameter set in the AutoWah_init function.
*/
void AutoWah_sweep(void) {
	double yout;
     
	if (!--counter) {
		if (!control) {
			bp_iir_setup(&H,(center_freq+=f_step));
			if (center_freq > max_freq) {
				control = 1;
			}
		}
		else if (control) {
			bp_iir_setup(&H,(center_freq-=f_step));
			if (center_freq == min_freq) {
				control = 0;
			}
		}

		counter = counter_limit; 
	}
}

/*************** AutoWah.h ****************************/

#ifndef __AUTOWAH_H__
#define __AUTOWAH_H__

extern void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step);
extern double AutoWah_process(int xin);
extern void AutoWah_sweep(void);

#endif

/************** Usage Example **************************/

#include "AutoWah.h"

void main(void) {
    short xin;
    short yout;
    AutoWah_init(2000,  /*Effect rate 2000*/
                 16000, /*Sampling Frequency*/
                 1000,  /*Maximum frequency*/
                 500,   /*Minimum frequency*/
                 4,     /*Q*/
                 0.707, /*Gain factor*/
                 10     /*Frequency increment*/
                 );   

    while(1) {
        if (new_sample_flag()) {
            /*When there's new sample at your ADC or CODEC input*/
            /*Read the sample*/
            xin = read_sample();
            /*Apply AutoWah_process function to the sample*/
            yout =AutoWah_process(xin);

            /*Send the output value to your DAC or codec output*/
            write_output(yout);
            /*Makes the LFO vary*/
            AutoWah_sweep();
        }                              
    }
}

Yet Another FIR design algorithm

Markus Nentwig ● August 21, 2011●3 comments ● Coded in Matlab

% ********************************************
% least-mean-squares FIR design algorithm
% Markus Nentwig, 2010-2011
% release 2011/8/22
% ********************************************
function LMSFIR()
    close all;
    
    h1 = demo1('basic'); compareDemo1WithRemez(h1);    
    % h1 = demo1('basicLMS'); disp('demo: convergence failure is on purpose to show LMS solution');
    % demo1('allpass');
    % demo1('stopband');
    % demo1('equalize');
    % demo1('nominalResponse');
    % demo1('multipassband');
    % demo1('complex');
    % demo1('rootRaisedCosineUpsampler');
    % demo1('componentModel');
    % demo1('componentModel2');
end

function h = demo1(nameOfDemo)
    dpar = struct();
    % parameters for a basic FIR lowpass filter design.
    % kept in a struct(), so that individual examples
    % can easily change them.    
    
    % sampling rate at the input of the filter
    dpar.inRate_Hz = 104e6;

    % number of physical FIR taps 
    % in a polyphase decimator, the number of internal
    % coefficients will be fDecim * fStages
    dpar.mStages = 36; 
    
    % end of passband. A single value will be internally
    % expanded to [-9e6, 9e6]. 
    % Asymmetric designs require
    % the complexValued = true option.
    % This 'default' passband can be omitted entirely, if passbands
    % are declared individually later
    dpar.req_passbandF_Hz = 9e6;
    
    % defines the maximum allowed ripple in the passband.
    dpar.req_passbandRipple_dB = 0.083;
    
    % as alternative to ripple, the in-band error
    % vector magnitude (EVM) can be limited
    % dpar.req_passbandEVM_dB = -44;
    
    % the passband specification may use multiple points
    % dpar.req_passbandF_Hz = [-1, 5e6, 6e6, 9e6];
    % dpar.req_passbandEVM_dB = [-44, -44, -34, -34];
    
    % start of default stopband. 
    % as with the passband, the default stopband can be omitted, 
    % if individual bands are placed later.
    dpar.req_stopbandF_Hz = 14e6;  
    dpar.req_stopbandMaxGain_dB = -30;

    % dpar.req_stopbandF_Hz = [14e6, 34e6];  
    % dpar.req_stopbandGainMax_dB = [-30, -20];
    
    % ********************************************
    % create a filter design object "design"
    % * access with LMSFIR_stage2 functions
    % * evaluate with LMSFIR_stage3 function    
    % ********************************************    

    switch nameOfDemo
      case 'basic'
        % ********************************************
        % simple filter using the parameters above
        % ********************************************
        design = LMSFIR_stage1_setup(dpar);
      
      case 'basicLMS'
        % ********************************************
        % LMS design for comparison:
        % Iterations are disabled
        % ********************************************
        dpar.nIter = 1;

        % balance in-band / out-of-band performance as needed
        dpar.inbandWeight_initValue = 5;
        dpar.outOfBandWeight_initValue = 1;
        
        design = LMSFIR_stage1_setup(dpar);
        
      case 'allpass'
        % ********************************************
        % allpass design Offset the nominal delay by 1/3
        % of a sample, compared to the "basic" example
        % (compare the impulse responses)
        % ********************************************
        dpar.delayOffset = 1/3; % signal arrives now earlier
        design = LMSFIR_stage1_setup(dpar);
        
      case 'stopband'        
        % ********************************************
        % Filter with added stopbands
        % ********************************************
        % the following features require more taps
        dpar.mStages = 48;
        
        % create filter design object
        design = LMSFIR_stage1_setup(dpar);
        
        % place a stopband from 14 to 16 MHz with -50 dB
        design = LMSFIR_stage2_placeStopband(...
            design, ...
            'f_Hz', [14e6, 16e6], ...
            'g_dB', [-50, -50]);
        
        % place another stopband from 16 to 28 MHz with 
        % -50 dB, linearly relaxing to -40 dB
        design = LMSFIR_stage2_placeStopband(...
            design, ...
            'f_Hz', [16e6, 28e6], ...
            'g_dB', [-50, -40]);

      case 'equalize'
        % ********************************************
        % Equalize the frequency response of another 
        % filter in the passband(s)
        % ********************************************
        % As an equalizer, this is rather inefficient with so much 
        % unused bandwidth. Should operate at the smallest possible BW instead.
        dpar.mStages = 52; 
        
        [ffilter_Hz, H] = getExampleLaplaceDomainFilter();

        % set the frequency points...
        dpar.equalizeFreq_Hz = ffilter_Hz;

        % ... and the filter response. The design routine will
        % use linear interpolation, therefore provide a sufficiently
        % dense grid.
        % Equalizing an asymmetric response requires 
        % the complexValued=true option, leading to a complex-valued
        % FIR filter.
        % The equalization function needs to be normalized. 
        % Otherwise, pass- and stopband targets will be offset 
        % by the gain mismatch.
        dpar.equalizeComplexGain = H;

        % as alternative to the complex gain, a magnitude response
        % can be given via an equalizeGain_dB argument.
        % dpar.equalizeGain_dB = 20*log10(abs(H));        

        % an asymmetric (non-linear-phase) impulse response may 
        % require a group delay that is not centered in the
        % FIR length.
        dpar.delayOffset = 2;
        
        design = LMSFIR_stage1_setup(dpar);
      case 'componentModel'
        % ********************************************
        % Create a FIR filter that approximates the passband behavior of 
        % the analog filter accurately, and gives a similar stopband rejection
        % 
        % The most straightforward way to model an infinite-impulse-response
        % lowpass is to simply sample the impulse response. However, it needs to be
        % cut to size (since the FIR filter has only finite length)
        % => Chopping it off destroys the out-of-band performance (=rectangular window)
        % => use a window function that trades off between passband accuracy and 
        %    stopband rejection
        % => or use the design example below.
        % ********************************************
        dpar.mStages = 52; 
        
        [ffilter_Hz, H] = getExampleLaplaceDomainFilter();

        % set the frequency points...
        dpar.nominalFreq_Hz = ffilter_Hz;
        dpar.nominalComplexGain = H;
        
        dpar.req_stopbandF_Hz = [15e6, 30e6, 55e6];
        dpar.req_stopbandMaxGain_dB = [-38, -80, -115];

        dpar.req_passbandF_Hz = 9e6;        

        % offset the impulse response, it is not centered
        dpar.delayOffset = 18;

        design = LMSFIR_stage1_setup(dpar);
      
      case 'componentModel2'
        % ********************************************
        % an extension of "componentModel1"
        % stopband rejection does not matter, but we need
        % phase-accurate modeling on a region of the stopband edge
        % ********************************************
        dpar.mStages = 80; % this won't be cheap...
      
        [ffilter_Hz, H] = getExampleLaplaceDomainFilter();
        dpar.nominalFreq_Hz = ffilter_Hz;
        dpar.nominalComplexGain = H;
        
        dpar.req_stopbandF_Hz = [ 16e6, 100e6];
        dpar.req_stopbandMaxGain_dB = [ -30, -30];
        
        dpar.req_passbandF_Hz = 9e6;        
        
        % offset the impulse response, it is not centered
        dpar.delayOffset = dpar.mStages / 2 - 8;
        
        design = LMSFIR_stage1_setup(dpar);
            
        % place a passband in the area on the slope that is to be modeled accurately
        design = LMSFIR_stage2_placePassband(...
            design, ...
            'f_Hz', [12e6, 16e6], ...
            'EVM_dB', [-40, -50] - 30); % nominal gain -40..-50 dB, -30 dBc EVM
        
      case 'nominalResponse'
        % ********************************************
        % Design a filter to a given frequency response
        % ********************************************    
        dpar.mStages = 50;
        % the frequency response is approximated in any
        % declared passband, but must be valid for any 
        % frequency to allow plotting.

        dpar.nominalFreq_Hz = [0, 3e6, 9e6, 1e9];
        dpar.nominalGain_dB = [0, 0, -6, -6]; 
        
        % instead, nominalComplexGain can be used
        % g = [0, 0, -3, -3]; 
        % dpar.nominalComplexGain = 10 .^ (g/20);
        design = LMSFIR_stage1_setup(dpar);
        
      case 'multipassband'
        % ********************************************
        % Design a filter with three passbands
        % ********************************************    
        dpar.mStages = 50;
        dpar = rmfield(dpar, 'req_passbandF_Hz');
        dpar = rmfield(dpar, 'req_passbandRipple_dB');
        
        design = LMSFIR_stage1_setup(dpar);

        design = LMSFIR_stage2_placePassband(...
            design, ...
            'f_Hz', [-2e6, 2e6], ...
            'EVM_dB', -45);
        
        design = LMSFIR_stage2_placePassband(...
            design, ...
            'f_Hz', [3e6, 7e6], ...
            'EVM_dB', [-45, -40]);

        design = LMSFIR_stage2_placeStopband(...
            design, ...
            'f_Hz', [11.8e6, 12.4e6], ...
            'g_dB', -70);

      case 'complex'
        % ********************************************
        % Design a complex-valued filter
        % ********************************************    
        % this is also an example for what can go wrong: 
        % In the unconstrained section around -40 MHz, the
        % frequency response is allowed to go berserk. Which
        % it does.
        % Solution: Place a "soft" stopband (for example at -5 dB)
        % in the "don't-care" regions and add a couple of taps.
        
        % remove passband from default parameters
        dpar = rmfield(dpar, 'req_passbandF_Hz');
        dpar = rmfield(dpar, 'req_passbandRipple_dB');

        % remove stopband from default parameters
        dpar = rmfield(dpar, 'req_stopbandF_Hz');
        dpar = rmfield(dpar, 'req_stopbandMaxGain_dB');
        
        dpar.complexValued = true;        
        design = LMSFIR_stage1_setup(dpar);

        design = LMSFIR_stage2_placeStopband(...
            design, ...
            'f_Hz', [-30e6, -16e6], ...
            'g_dB', -50);

        design = LMSFIR_stage2_placePassband(...
            design, ...
            'f_Hz', [-8e6, -2e6], ...
            'EVM_dB', -45);
        
        design = LMSFIR_stage2_placeStopband(...
            design, ...
            'f_Hz', [3e6, 40e6], ...
            'g_dB', [-30, -50]);

      case 'rootRaisedCosineUpsampler'
        % ********************************************
        % root-raised cosine upsampling filter for WCDMA transmitter
        % The input chip stream arrives at 3.84 Msps, using the 
        % full bandwidth. 
        % Before the filter, it is upsampled (zero insertion) to 
        % 7.68 Msps. 
        % The filter applies RRC-filtering with 1.22 rolloff.
        % ********************************************
        % calculate nominal RRC response for lookup table / linear
        % interpolation
        f_Hz = logspace(1, 8, 10000); f_Hz(1) = -1;
        c = abs(f_Hz * 2 / 3.84e6);
        c = (c-1)/(0.22); % -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);
        
        % ********************************************
        % once the targets are achieved, use the remaining
        % 'degrees of freedom' for least-squares optimization.
        % The LMS solver will improve, where it is 'cheapest'
        % The parameters are not a real-world design 
        % (0.5 percent EVM => -46 dB)
        % ********************************************    
        ci = 0;
        % ci = 10; % for comparison: force equiripple 
        
        dpar = struct(...
            'inRate_Hz', 3.84e6, ...
            'fInterp', 2, ...
            'mStages', 45, ...
            'req_passbandF_Hz', 3.84e6 * 1.22 / 2, ...
            'req_passbandEVM_dB', -46, ...
            'req_stopbandF_Hz', 2.46e6, ...
            'req_stopbandMaxGain_dB', -50, ...
            'nominalFreq_Hz', f_Hz, ...
            'nominalGain_dB', 20*log10(RRC_h + 1e-19), ...
            'convergedIterations', ci);
        
        design = LMSFIR_stage1_setup(dpar);
        [h, status] = LMSFIR_stage3_run(design); 
        % save('hRRC_upsampler.txt', 'h', '-ascii');
        disp(status);        
        
      otherwise
        assert(false);
    end % switch nameOfDemo
    
    % ********************************************
    % Design the filter
    % ********************************************    
    % h is the impulse response (FIR tap coefficients).
    [h, status] = LMSFIR_stage3_run(design); 
    disp(status);    
end

function compareDemo1WithRemez(hLMS)
% identical target settings to demo1 "basic".
% note, the demo uses targets that are exactly "on the edge"
% what the algorithm can achieve. This results in an equiripple-
% design that can be compared with remez().
% If the targets are too loosely set, pass- and stopband quality
% start to "sag" in the band middle (LMS solution => lowest overall
% error, the optimizer improves where it's "cheapest").
    
    r_Hz = 104e6;
    m = 35; % definition differs by 1
    f = [0 9e6 14e6 r_Hz/2] / (r_Hz/2);
    a = [1 1 0 0];
    ripple_dB = 0.1;
    att_dB = 30;
    err1 = 1 - 10 ^ (-ripple_dB / 20);
    err2 = 10 ^ (-att_dB / 20);
    w = [1/err1 1/err2];

    % get remez design impulse response
    hRemez = remez(m, f, a, w);
    
    figure(); hold on;
    handle = plot(hLMS, 'b+'); set(handle, 'lineWidth', 3);
    plot(hRemez, 'k+'); set(handle, 'lineWidth', 2);
    legend('this algorithm', 'Remez');
    title('comparison with Remez design (optimum equiripple)');
end

% Gets the frequency response of an "analog" (Laplace-domain) filter.
% => Chebyshev response
% => 6th order
% => 1.2 dB ripple
% => cutoff frequency at 10 MHz
% returns
% f_Hz: list of frequencies
% H: complex-valued H(f_Hz)
function [f_Hz, H] = getExampleLaplaceDomainFilter()
    [IIR_b, IIR_a] = cheby1(6, 1.2, 1, 's');
    
    % evaluate it on a wide enough frequency range
    f_Hz = logspace(1, 10, 1000); f_Hz(1) = -1;
    
    % Laplace domain operator for normalized frequency
    fCutoff_Hz = 10e6;
    s = 1i * f_Hz / fCutoff_Hz;
    
    % polynomial in s
    H_num = polyval(IIR_b, s);
    H_denom = polyval(IIR_a, s);
    H = H_num ./ H_denom;
end

% === LMSFIR_xyz "API" functions ===

% ********************************************
% LMSFIR_stagex_... functions to interact with design 'object'
% to be executed in 'stage'-order
% ********************************************
function d = LMSFIR_stage1_setup(varargin)
    p = varargin2struct(varargin);
    d = struct();

    % number of frequency points. Increase to improve accuracy.
    % Frequencies are quantized to +/- rate / (2 * nSamples) 
    d.nSamples = 1024;
    
    % default polyphase interpolation: none 
    d.fInterp = 1;
    
    % default polyphase decimation: none
    d.fDecim = 1;

    % max. number of iterations
    d.nIter = 100; 

    % for pure LMS solution, set nIter to 1 and change the weights below as needed
    d.inbandWeight_initValue = 1;
    d.outOfBandWeight_initValue = 1;
    
    % abort when the iteration weights grow too large. 
    % This happens when targets are impossible.
    % The result may still be meaningful, though.
    d.abortWeight = 1e12;

    % keep iterating, if the targets are reached. 
    % Once the "equi"-ripple iteration has brought all peaks to an acceptable level, 
    % the LMS solver will use the remaining "degrees of freedom" for a LMS optimization.
    % The solver improves "where it is easy / cheap". This results in sloped
    % stopbands and "drooping" EVM in passbands.
    % Often, LMS is the best choice => set converged iterations to 0.
    d.convergedIterations = 10;

    % for a complex-valued filter, use "true".
    % With a real-valued design, user input is only evaluated for positive frequencies!
    d.complexValued = false;
    
    % by default, the basis waveforms given to the optimizer are
    % within a delay range of +/- half the FIR length.
    % For nonlinear phase types (equalization / nominal frequency 
    % response), this may be suboptimal.
    % Meaningful values shouldn't exceed +/- half the number of
    % coefficients in the impulse response.
    d.delayOffset = 0;

    % copy parameters
    fn = fieldnames(p);
    for ix = 1:size(fn, 1)
        key = fn{ix};
        d.(key) = p.(key);           
    end

    % frequency base over FFT range
    fb = 0:(d.nSamples - 1);
    fb = fb + floor(d.nSamples / 2);
    fb = mod(fb, d.nSamples);
    fb = fb - floor(d.nSamples / 2);
    fb = fb / d.nSamples; % now [0..0.5[, [-0.5..0[
    fb = fb * d.inRate_Hz * d.fInterp;
    d.fb = fb;
    
    % in real-valued mode, negative frequencies are treated as 
    % positive, when 'user input' is evaluated
    if d.complexValued
        d.fbUser = fb;
    else
        d.fbUser = abs(fb);
    end
    
    % ********************************************
    % target settings. Those will be modified by 
    % LMSFIR_stage2_xyz()
    % ********************************************
    % initial value of NaN indicates: all entries are unset
    d.errorSpecBinVal_inband_dB = zeros(size(d.fb)) + NaN;
    d.errorSpecBinVal_outOfBand_dB = zeros(size(d.fb)) + NaN;
    
    % ********************************************
    % process req_passband requirement
    % needs to be done at stage 1, because it is 
    % used for 'gating' with the tightenExisting / 
    % relaxExisting options
    % ********************************************
    if isfield(d, 'req_passbandF_Hz')  

        par = struct('onOverlap', 'error');
        if isfield(d, 'req_passbandRipple_dB')
            par.ripple_dB = d.req_passbandRipple_dB;
        end
        if isfield(d, 'req_passbandEVM_dB')
            par.EVM_dB = d.req_passbandEVM_dB;
        end
        par.f_Hz = d.req_passbandF_Hz;

        d = LMSFIR_stage2_placePassband(d, par);
    end % if req_passbandF_Hz

    % ********************************************
    % process req_stopband requirement
    % needs to be done at stage 1, because it is 
    % used for 'gating' with the tightenExisting / 
    % relaxExisting options
    % ********************************************
    if isfield(d, 'req_stopbandF_Hz')  
        f_Hz = d.req_stopbandF_Hz;
        g_dB = d.req_stopbandMaxGain_dB;
        
        % extend to infinity
        if isscalar(f_Hz)
            f_Hz = [f_Hz 9e19]; 
            g_dB = [g_dB g_dB(end)];
        end
        
        d = placeBand...
            (d, ...
             'f_Hz', f_Hz, 'g_dB', g_dB, ...
             'type', 'stopband', ...
             'onOverlap', 'tighten');
    end    
    
    % ********************************************
    % plot management
    % ********************************************
    d.nextPlotIx = 700;
end

function d = LMSFIR_stage2_placeStopband(d, varargin)
    p = varargin2struct(varargin);
    
    % shorthand notation g_dB = -30; f_Hz = 9e6;
    % extend fixed passband to positive infinity
    if isscalar(p.f_Hz)
        assert(p.f_Hz > 0);
        p.f_Hz = [p.f_Hz 9e99];
    end

    if isscalar(p.g_dB)
        p.g_dB = ones(size(p.f_Hz)) * p.g_dB;
    end
    
    % default action is to use the stricter requirement
    if ~isfield(p, 'onOverlap')
        p.onOverlap = 'tighten';
    end        
    
    d = placeBand(d, 'type', 'stopband', p);
end

function d = LMSFIR_stage2_placePassband(d, varargin)
    p = varargin2struct(varargin);
    
    % default action is to use the stricter requirement
    if ~isfield(p, 'onOverlap')
        p.onOverlap = 'tighten';
    end        
    
    % translate ripple spec to error
    if isfield(p, 'ripple_dB')
        assert(p.ripple_dB > 0);
        eSamplescale = 10 ^ (p.ripple_dB / 20) - 1;            
        EVM_dB = 20*log10(eSamplescale);
    end
    
    if isfield(p, 'EVM_dB')
        EVM_dB = p.EVM_dB;
    end

    % convert scalar to two-element vector
    if isscalar(EVM_dB)
        EVM_dB = [EVM_dB EVM_dB];
    end
    
    % *** handle f_Hz ***
    f_Hz = p.f_Hz;
    
    % extend to 0 Hz
    if isscalar(f_Hz)
        f_Hz = [0 f_Hz];
    end
    
    % *** create the passband ***
    d = placeBand(d, ...
                  'type', 'passband', ...
                  'f_Hz', f_Hz, ...
                  'g_dB', EVM_dB, ...
                  'onOverlap', p.onOverlap);
end

% ********************************************
% the filter design algorithm
% h: impulse response
% status: converged or not
% note that even if convergence was not reached, 
% the resulting impulse response is "the best we
% can do" and often meaningful.
% ********************************************
function [h, status] = LMSFIR_stage3_run(d)
    1;    
    
    % mTaps is number of physical FIR stages
    % m is number of polyphase coefficients
    d.m = d.mStages * d.fInterp;
    
    % masks flagging pass-/stopband frequencies
    mask_inband = NaN_to_0_else_1(d.errorSpecBinVal_inband_dB);
    mask_outOfBand = NaN_to_0_else_1(d.errorSpecBinVal_outOfBand_dB);
    
    % sanity check... (orthogonality of wanted and unwanted component)
    assert(sum(mask_inband) > 0, 'passband is empty');
    assert(sum(mask_inband .* mask_outOfBand) == 0, ...
           'passband and stopband overlap');
    
    % ********************************************
    % start with flat passband signals at input and output of filter
    % those will become the input to the LMS solver.
    % ********************************************
    sigSolverAtInput_fd = mask_inband;
    sigSolverAtOutput_fd = sigSolverAtInput_fd;

    % ********************************************
    % for even-sized FFT length, there is one bin at the
    % Nyquist limit that gives a [-1, 1, -1, 1] time domain
    % waveform. It has no counterpart with opposite frequency
    % sign and is therefore problematic (time domain component
    % cannot be delayed). 
    % Don't assign any input power here.
    % ********************************************
    if mod(d.nSamples, 2) == 0
        ixNyquistBin = floor(d.nSamples/2) + 1;
        sigSolverAtInput_fd(ixNyquistBin) = 0;
        sigSolverAtOutput_fd(ixNyquistBin) = 0;
    end
    
    if isfield(d, 'equalizeFreq_Hz')
        % ********************************************
        % Filter equalizes a given passband frequency response
        % ********************************************
        if isfield(d, 'equalizeGain_dB')
            cgain = 10 .^ (equalizeGain_dB / 20);
        else
            cgain = d.equalizeComplexGain;
        end
        d.Heq = evaluateFilter(d.fb, d.equalizeFreq_Hz, cgain, d.complexValued);
        assert(isempty(find(isnan(d.Heq), 1)), ...
               ['equalizer frequency response interpolation failed. ' ...
                'Please provide full range data for plotting, even if it does not ', ...
                'affect the design']);
        
        % ********************************************
        % apply frequency response to input signal. 
        % The LMS solver will invert this response
        % ********************************************
        sigSolverAtInput_fd = sigSolverAtInput_fd .* d.Heq;
    end

    if isfield(d, 'nominalFreq_Hz')
        % ********************************************
        % (equalized) filter matches a given passband frequency response
        % ********************************************
        if isfield(d, 'nominalGain_dB')
            cgain = 10 .^ (d.nominalGain_dB / 20);
        else
            cgain = d.nominalComplexGain;
        end
        
        d.Hnom = evaluateFilter(d.fb, d.nominalFreq_Hz, cgain, d.complexValued);
        assert(isempty(find(isnan(d.Hnom), 1)), ...
               ['nominal frequency response interpolation failed. ' ...
                'Please provide full range data for plotting, even if it does not ', ...
                'affect the design']);
        
        % ********************************************
        % apply frequency response to output signal. 
        % The LMS solver will adapt this response
        % ********************************************
        sigSolverAtOutput_fd = sigSolverAtOutput_fd .* d.Hnom;
    end
    
    % ********************************************
    % compensate constant group delay from equalizer and nominal
    % frequency response. This isn't optimal, but it is usually
    % a good starting point (use delayOffset parameter)
    % ********************************************
    [coeff, ref_shiftedAndScaled, deltaN] = fitSignal_FFT(...
        ifft(sigSolverAtInput_fd), ifft(sigSolverAtOutput_fd));

    % the above function also scales for best fit. This is not desired here, instead
    % let the LMS solver match the gain. Simply scale it back:
    ref_shifted = ref_shiftedAndScaled / coeff;
    sigSolverAtOutput_fd = fft(ref_shifted);
    
    if false
        % ********************************************
        % plot time domain waveforms (debug)
        % ********************************************
        figure(76); hold on;
        plot(fftshift(abs(ifft(sigSolverAtOutput_fd))), 'k');
        plot(fftshift(abs(ifft(sigSolverAtInput_fd))), 'b');
        title('time domain signals');
        legend('reference (shifted)', 'input signal');
    end
    
    % ********************************************
    % main loop of the design algorithm
    % => initialize weights
    % => loop
    % =>   design optimum LMS filter that transforms weighted input 
    %      into weighted output
    % =>   adapt weights 
    % => iterate
    % ********************************************
    % at this stage, the input to the algorithm is as follows:
    % => errorSpec for in-band and out-of-band frequencies
    %    (masks are redundant, can be derived from above)
    % => LMS_in_fd and
    % => LMS_out_fd: Signals that are given to the LMS solver.
    %   Its task is: "design a FIR filter that transforms LMS_in_fd into LMS_out_fd".
    
    % initialize weights
    outOfBandWeight = mask_outOfBand * d.outOfBandWeight_initValue;
    inbandWeight = mask_inband * d.inbandWeight_initValue;
    
    status = '? invalid ?';
    hConv = [];
    remConvIter = d.convergedIterations;
    for iter=1:d.nIter
        % inband weight is applied equally to both sides to shape the error
        % out-of-band weight is applied to the unwanted signal
        LMS_in_fd = sigSolverAtInput_fd .* inbandWeight...
            + mask_outOfBand .* outOfBandWeight; 
        
        LMS_out_fd = sigSolverAtOutput_fd .* inbandWeight;
        
        % ********************************************
        % cyclic time domain waveforms from complex spectrum
        % ********************************************
        LMS_in_td = ifft(LMS_in_fd); 
        LMS_out_td = ifft(LMS_out_fd);
        
        % ********************************************
        % construct FIR basis (output per coeffient)
        % time domain waveforms, shifted according to each FIR tap
        % ********************************************
        basis = zeros(d.m, d.nSamples);

        % introduce group delay target
        ix1 = -d.m/2+0.5 + d.delayOffset;
        ix2 = ix1 + d.m - 1;
        rowIx = 1;
        for ix = ix1:ix2 % index 1 appears at ix1                         
            basis(rowIx, :) = FFT_delay(LMS_in_td, ix);
            rowIx = rowIx + 1;
        end
        
        if d.complexValued
            rightHandSide_td = LMS_out_td;
        else
            % use real part only
            basis = [real(basis)];
            rightHandSide_td = [real(LMS_out_td)]; 

            pRp = real(rightHandSide_td) * real(rightHandSide_td)' + eps;
            pIp = imag(rightHandSide_td) * imag(rightHandSide_td)';
            assert(pIp / pRp < 1e-16, ...
                   ['got an imaginary part where there should be none. ', ...
                    'uncomment the following lines, if needed']);
            % if designing a real-valued equalizer for a complex-valued frequency response,
            % use the following to solve LMS over the average:
            % basis = [real(basis) imag(basis)];
            % rightHandSide_td = [real(LMS_out_td), imag(LMS_out_td)]; 
        end
        
        % ********************************************
        % LMS solver
        % find a set of coefficients that scale the 
        % waveforms in "basis", so that their sum matches
        % "rightHandSide_td" LMS-optimally
        % ********************************************
        pbasis = pinv(basis .');
        h = transpose(pbasis * rightHandSide_td .');
        
        % pad impulse response to n
        irIter = [h, zeros(1, d.nSamples-d.m)];

        % undo the nominal group delay
        irIter = FFT_delay(irIter, ix1);
        HIter = fft(irIter);
        
        % ********************************************
        % filter test signal
        % ********************************************
        eq_fd = sigSolverAtInput_fd .* HIter; 

        % ********************************************
        % subtract actual output from targeted output
        % results in error spectrum
        % ********************************************
        err_fd = sigSolverAtOutput_fd - eq_fd;
        err_fd = err_fd .* mask_inband; % only in-band matters
        EVM_dB = 20*log10(abs(err_fd)+1e-15);
        
        % ********************************************
        % out-of-band leakage
        % ********************************************
        leakage_dB = 20*log10(abs(HIter .* mask_outOfBand + 1e-15));

        % ********************************************
        % compare achieved and targeted performance 
        % ********************************************
        deltaLeakage_dB = leakage_dB - d.errorSpecBinVal_outOfBand_dB;        
        deltaEVM_dB = EVM_dB - d.errorSpecBinVal_inband_dB;
        
        % ********************************************
        % find bins where performance should be improved
        % or relaxed
        % ********************************************
        ixImprLeakage = find(deltaLeakage_dB > 0);
        ixImprEVM = find(deltaEVM_dB > 0);
        ixRelLeakage = find(deltaLeakage_dB < -3);
        ixRelEVM = find(deltaEVM_dB < -3);        
        
        status = 'iteration limit reached';
        if isempty(ixImprLeakage) && isempty(ixImprEVM)
            % both targets met. Convergence!
            if remConvIter > 0
                remConvIter = remConvIter - 1;
                status = 'converged once, now trying to improve';
                hConv = h;
            else
                status = 'converged';
                break;
            end
        end
        
        % ********************************************
        % improve / relax in-band and out-of-band
        % ********************************************
        if ~isempty(ixImprLeakage)
            % tighten out-of-band
            outOfBandWeight(ixImprLeakage) = outOfBandWeight(ixImprLeakage)...
                .* 10 .^ ((deltaLeakage_dB(ixImprLeakage) + 0.1) / 20);
        end       
        
        if ~isempty(ixRelLeakage)
            % relax out-of-band
            outOfBandWeight(ixRelLeakage) = outOfBandWeight(ixRelLeakage)...
                .* 10 .^ (deltaLeakage_dB(ixRelLeakage) / 3 / 20); 
        end
        
        if ~isempty(ixImprEVM)
            % tighten in-band
            inbandWeight(ixImprEVM) = inbandWeight(ixImprEVM)...
                .* 10 .^ ((deltaEVM_dB(ixImprEVM) + 0.01) / 20); 
        end
        if ~isempty(ixRelEVM)
            % relax in-band
            inbandWeight(ixRelEVM) = inbandWeight(ixRelEVM)...
                .* 10 .^ (deltaEVM_dB(ixRelEVM) / 2 / 20); 
        end
        
        if max([inbandWeight, outOfBandWeight] > d.abortWeight)
            status = 'weight vector is diverging';
            break;
        end
    end % for iter
    
    % ********************************************
    % recover from convergence failure after convergence
    % during improvement phase
    % ********************************************
    if ~strcmp(status, 'converged')
        if ~isempty(hConv)
            h = hConv;
            status = 'converged';
        end
    end
    
    if true
        % ********************************************
        % plot impulse response
        % ********************************************
        if d.complexValued
            figure(); hold on; 
            stem(real(h), 'k');
            stem(imag(h), 'b');
            legend('real(h)', 'imag(h)');
        else
            figure(); hold on; 
            stem(h);
            legend('h');
        end
        title('impulse response');
    end
    
    
    % ********************************************
    % plot frequency response
    % ********************************************
    inbandBins = find(mask_inband);
    outOfBandBins = find(mask_outOfBand);
    d=doPlotStart(d, ['Frequency response (Status:', status, ')']);
    d=doPlotH(d, HIter, 'b', '|H_{design}(f)|', 2);
    
    handle = plot(d.fb(outOfBandBins), d.errorSpecBinVal_outOfBand_dB(outOfBandBins), 'b+');
    set(handle, 'markersize', 2);
    d=addLegend(d, 'req. stopband');

    d = doPlot_dB(d, EVM_dB, 'r', 'error');
    
    handle = plot(d.fb(inbandBins), d.errorSpecBinVal_inband_dB(inbandBins), 'r+');
    set(handle, 'markersize', 2);
    d=addLegend(d, 'req. passband error');
    
    d=doPlotEnd(d);
    ylim([-100, 10]);
    
    if false
        % ********************************************
        % plot constraint signal and weights
        % ********************************************
        figure(31); grid on; hold on;
        handle = plot(fftshift(d.fb), fftshift(20*log10(mask_outOfBand))); 
        set(handle, 'lineWidth', 3);
        
        x = d.fb; y = 20*log10(inbandWeight / max(inbandWeight));
        handle = plot(x(inbandBins), y(inbandBins), 'k+'); set(handle, 'lineWidth', 3);
        
        x = d.fb;
        y = 20*log10(outOfBandWeight / max(outOfBandWeight));
        handle = plot(x(outOfBandBins), y(outOfBandBins), 'b+'); set(handle, 'lineWidth', 3); 
        
        xlabel('f/Hz'); ylabel('dB');
        ylim([-80, 40]);
        legend('constraint signal', 'in-band weight', 'out-of-band weight');
        title('weighting factor (normalized to 0 dB)');
    end
    
    hasEq = isfield(d, 'Heq');
    hasNom = isfield(d, 'Hnom');
    if hasEq || hasNom
        % ********************************************
        % plot equalization / nominal target
        % ********************************************
        d=doPlotStart(d, 'equalization / nominal target');
        d=doPlotH(d, HIter, 'b', '|H_{design}(f)|', 2);
        if hasEq
            d=doPlotH(d, d.Heq, 'k', '|H_{eq}(f)| to equalize (invert)');
            eqR = HIter .* d.Heq;
            d=doPlotH(d, eqR, 'c', '|H_{design}(f)H_{eq}(f)|', 2);
            handle = plot(d.fb(inbandBins), ...
                          20*log10(abs(eqR(inbandBins)) + 1e-15), 'c*');
            set(handle, 'markersize', 3);
            d=addLegend(d, '|H_{eq}(in-band)');
        end
        if hasNom
            d = doPlotH(d, d.Hnom, 'g', '|H_{nom}|', 2);
            handle = plot(d.fb(inbandBins), ...
                          20*log10(abs(HIter(inbandBins)) + 1e-15), 'b*');
            set(handle, 'markersize', 3);
            d=addLegend(d, '|H_{design}(f)H_{eq}(f) in-band');
        end
        d=doPlotEnd(d);
        % set y-range
        ymax = 20*log10(max(abs(HIter)));
        ylim([-50, ymax+3]);
    end
end

% === LMSFIR helper functions ===

% evaluates frequency response f_dB; g_Hz at fb
% the return value will contain NaN for out-of-range entries
% in fb
function binVal = buildBinVal(varargin)
    p = varargin2struct(varargin);
    
    f_Hz = p.f_Hz;
    g_dB = p.g_dB;

    % shorthand notation f = [f1, f2]; g = -30;
    if isscalar(g_dB)
        g_dB = ones(size(f_Hz)) * g_dB;
    end
    
    % tolerate sloppy two-argument definition
    if size(f_Hz, 2) == 2 && f_Hz(1) > f_Hz(2)
        f_Hz = fliplr(f_Hz);
        g_dB = fliplr(g_dB);
    end
    
    binVal = interp1(f_Hz, g_dB, p.fbUser, 'linear');
end

function d = placeBand(d, varargin)
    p = varargin2struct(varargin);
    
    % create requirements vector
    binVal = buildBinVal('f_Hz', p.f_Hz, ...
                         'g_dB', p.g_dB, ...
                         'fbUser', d.fbUser);
    
    % look up requirements vector from design object
    switch p.type
      case 'passband'
        fn = 'errorSpecBinVal_inband_dB';
      case 'stopband'
        fn = 'errorSpecBinVal_outOfBand_dB';
      otherwise
        assert(false);
    end
    designObject_binVal = d.(fn);

    % check overlap
    if strcmp(p.onOverlap, 'error')
        m1 = NaN_to_0_else_1(designObject_binVal);
        m2 = NaN_to_0_else_1(binVal);
        assert(isempty(find(m1 .* m2, 1)), ...
               ['newly declared band overlaps existing band, '...
                'which was explicitly forbidden by onOverlap=error']);
        p.onOverlap = 'tighten'; % there won't be overlap,
                                 % merging is dummy operation
    end
    
    % merging rules
    switch p.onOverlap
      case 'tighten'
        logicOp = 'or';
        valueOp = 'min';
      case 'relax'
        logicOp = 'or';
        valueOp = 'max';
      case 'tightenExisting'
        logicOp = 'and';
        valueOp = 'min';
      case 'relaxExisting'
        logicOp = 'and';
        valueOp = 'max';
      otherwise
        assert(false);
    end
    
    % merge requirements tables
    binValMerged = mergeBinVal(...
        'binVal1', designObject_binVal, ...
        'binVal2', binVal, ...
        'logicalOperator', logicOp, ...
        'valueOperator', valueOp);

    % assign new requirements table
    d.(fn) = binValMerged;    
end

function r = NaN_to_0_else_1(vec)
    r = zeros(size(vec));
    % logical indexing, instead of r(find(~isnan(vec))) = 1;
    r(~isnan(vec)) = 1; 
end

function binVal = mergeBinVal(varargin)
    p = varargin2struct(varargin);
    
    % region where first argument is defined
    mask1 = NaN_to_0_else_1(p.binVal1);

    % region where second argument is defined    
    mask2 = NaN_to_0_else_1(p.binVal2);

    % region where result will be defined    
    switch(p.logicalOperator)
      case 'or'
        mask = mask1 + mask2;
      case 'and'
        mask = mask1 .* mask2;
      otherwise
        assert(false);
    end 
    ix = find(mask);

    % merge into result
    binVal = zeros(size(p.binVal1)) + NaN;
    switch(p.valueOperator)
      case 'min'
        % note: The function min/max ignore NaNs (see "min" man 
        % page in Matlab)
        % if one entry is NaN, the other entry will be returned
        binVal(ix) = min(p.binVal1(ix), p.binVal2(ix));
      case 'max'
        binVal(ix) = max(p.binVal1(ix), p.binVal2(ix));
      otherwise
        assert(false);
    end
end        

% evaluates [f / gain] filter specification on the frequency grid
function H = evaluateFilter(f_eval, modelF, modelH, complexValued)

    oneSided = false;
    if ~complexValued
        oneSided = true;
    else
        if min(modelF) > min(f_eval)
            disp(['Warning: Filter model does not contain (enough) negative frequencies. ', ...
                  'assuming symmetric H(f) / real-valued h(t)']);
            oneSided = true;
        end
    end

    if oneSided
        f_evalOrig = f_eval;
        f_eval = abs(f_eval);
    end
    
    H = interp1(modelF, modelH, f_eval, 'linear');

    if oneSided
        % enforce symmetry (=> real-valued impulse response)
        logicalIndex = (f_evalOrig < 0);
        H(logicalIndex) = conj(H(logicalIndex));
    end
end

function [d, handle] = doPlotH(d, H, spec, legEntry, linewidth)
    handle = plot(fftshift(d.fb), fftshift(20*log10(abs(H)+1e-15)), spec);
    d = addLegend(d, legEntry);
    if exist('linewidth', 'var')
        set(handle, 'lineWidth', linewidth);
    end
end

function [d, handle] = doPlot_dB(d, H, spec, legEntry, linewidth)
    handle = plot(fftshift(d.fb), fftshift(H), spec);
    d.legList{size(d.legList, 2) + 1} = legEntry;
    if exist('linewidth', 'var')
        set(handle, 'lineWidth', linewidth);
    end
end

function d = doPlotStart(d, plotTitle)
    figure(d.nextPlotIx); 
    title(plotTitle);
    grid on; hold on;
    d.nextPlotIx = d.nextPlotIx + 1;
    d.legList = {};
end

function d = doPlotEnd(d)
    legend(d.legList);
    xlabel('f/Hz');
    ylabel('dB');
end

function d = addLegend(d, legEntry)
    d.legList{size(d.legList, 2) + 1} = legEntry;
end

% === general-purpose library functions ===

% handling of function arguments
% someFun('one', 1, 'two', 2, 'three', 3) => struct('one', 1, 'two', 2, 'three', 3)
% a struct() may appear in place of a key ('one') and gets merged into the output.
function r = varargin2struct(arg)
    assert(iscell(arg));
    
    switch(size(arg, 2))
      case 0 % varargin was empty
        r=struct();    
      case 1 % single argument, wrapped by varargin into a cell list 
        r=arg{1}; % unwrap
        assert(isstruct(r));
      otherwise
        r=struct();
        % iterate through cell elements
        ix=1;
        ixMax=size(arg, 2);
        while ix <= ixMax
            e=arg{ix};
            if ischar(e)
                % string => key/value. The next field is a value
                ix = ix + 1;
                v = arg{ix};
                r.(e) = v;
            elseif isstruct(e)
                names = fieldnames(e);
                assert(size(names, 2)==1); % column
                for ix2 = 1:size(names, 1) 
                    k = names{ix2};
                    v = e.(k);
                    r.(k) = v;
                end
            else 
                disp('invalid token in vararg handling. Expecting key or struct. Got:');
                disp(e);
                assert(false)
            end
            ix=ix+1;
        end % while
    end % switch
end

function sig = FFT_delay(sig, nDelay)
    sig = fft(sig); % to frequency domain
    nSigSamples = size(sig, 2);
    binFreq=(mod(((0:nSigSamples-1)+floor(nSigSamples/2)), nSigSamples)-floor(nSigSamples/2));
    phase = -2*pi*nDelay / nSigSamples .* binFreq;
    rot = exp(1i*phase);
    if mod(nSigSamples, 2)==0
        % even length - bin at Nyquist limit
        rot(nSigSamples/2+1)=cos(phase(nSigSamples/2+1));
    end
    sig = sig .* rot;
    sig = ifft(sig); % to time domain
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

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;

Remez (FIR design) weights from requirements

Markus Nentwig ● August 19, 2011●1 comment ● Coded in Matlab

% *********************************************
% Weighting factors for Remez equiripple FIR design
% M. Nentwig
% *********************************************
close all; clear all;

f = [0:9]/10; % normalized frequency 0..1

% a = nominal |H(f)| on a linear scale (sample amplitude)
% 1 1     : passband from frequency 0 to 0.1 
% 0 0     : stopband from frequency 0.2 to 0.3 
% 0.5 0.5 : passband with -6 dB (0.5 amplitude) from frequency 0.4 to 0.5 
% 0 0     : stopband from frequency 0.6 to 0.7
% 1 1     : passband from frequency 0.8 to 0.9
% other frequency ranges are "don't care" areas.
a = [1 1 0 0 0.5 0.5 0 0 1 1];

% *********************************************
% design specification for each band
% *********************************************
ripple1_dB = 0.3;
att2_dB = 60;
ripple3_dB = 0.2;
att4_dB = 70;
ripple5_dB = 0.1;

% *********************************************
% error for each band, on a linear scale
% *********************************************
% note: passband 3 has a nominal gain of 0.5 in 'a'.
% the allowed error of |H(f)| scales accordingly.
err1 = 1 - 10 ^ (-ripple1_dB / 20);
err2 = 10 ^ (-att2_dB / 20);
err3 = (1 - 10 ^ (-ripple3_dB / 20)) * 0.5; 
err4 = 10 ^ (-att4_dB / 20);
err5 = 1 - 10 ^ (-ripple5_dB / 20);

% the weight of each band is chosen as the inverse of the targeted error
% (stricter design target => higher weight).
% we could normalize each entry in w to (err1+err2+err3+err4+err5)
% but it would appear as common factor in all entries and therefore make no difference.
w = [1/err1 1/err2 1/err3 1/err4 1/err5];

% filter order
% this design problem is 'tweaked' so that the resulting filter is (quite) exactly on target
% for the given n, which can be changed only in steps of 1.
% Note that increasing order by 1 can make the filter worse due to even / odd 
% number of points (impulse response symmetry)
n = 52;

% *********************************************
% Run Remez / Parks McClellan filter design
% *********************************************
h = remez(n, f, a, w);
% force row vector (OctaveForge and Matlab's "remez" differ)
h = reshape(h, 1, prod(size(h))); 
% *********************************************
% Plot the results
% *********************************************
figure(1); hold on;

% zero-pad the impulse response to set frequency resolution
% of the FFT
h = [h, zeros(1, 1000)];

% frequency base
npts = size(h, 2);
fbase = (0:npts-1)/npts; 

plot(fbase, 20*log10(abs((fft(h)+1e-15))), 'b');
xlim([0, 0.5]);

% *********************************************
% sketch the requirements
% *********************************************
e = [1 1];
plot([f(1), f(2)]/2, e * ripple1_dB, 'k');
plot([f(1), f(2)]/2, e * -ripple1_dB, 'k');
plot([f(3), f(4)]/2, e * -att2_dB, 'k');
plot([f(5), f(6)]/2, e * ripple3_dB - 6.02, 'k');
plot([f(5), f(6)]/2, e * -ripple3_dB - 6.02, 'k');
plot([f(7), f(8)]/2, e * -att4_dB, 'k');
plot([f(9), f(10)]/2, e * ripple5_dB, 'k');
plot([f(9), f(10)]/2, e * -ripple5_dB, 'k');
xlabel('normalized frequency 0..1');
ylabel('dB');

Previous 21 222324 25 26 Next

vectors_intr.asm - External Interruption Configuration

David Valencia ● February 2, 2011 ● Coded in ASM for the TI C67x

*Vectors_intr.asm Vector file for interrupt INT11 
   .global _vectors     ;global symbols 
   .global _c_int00 
   .global _vector1 
   .global _vector2 
   .global _vector3 
   .global _c_int04     ; símbolo para EXT_INT4 
   .global _vector5 
   .global _vector6 
   .global _vector7 
   .global _vector8 
   .global _vector9    
   .global _vector10  
   .global _c_int11        ;for INT11 
   .global _vector12   
   .global _vector13    
   .global _vector14 
   .global _vector15 
 
   .ref _c_int00        ;entry address 
 
VEC_ENTRY .macro addr      ;macro for ISR 
    STW   B0,*--B15 
    MVKL  addr,B0 
    MVKH  addr,B0 
    B     B0 
    LDW   *B15++,B0 
    NOP   2 
    NOP    
    NOP    
   .endm 
 
_vec_nmi: 
  B    NRP 
  NOP  5 
 
_vec_dummy: 
  B    IRP 
  NOP  5 
 
 .sect ".vecs"        ;aligned IST section 
 .align 1024 
_vectors: 
_vector0:   VEC_ENTRY _c_int00     ;RESET 
_vector1:   VEC_ENTRY _vec_nmi   ;NMI 
_vector2:   VEC_ENTRY _vec_dummy    ;RSVD 
_vector3:   VEC_ENTRY _vec_dummy 
_vector4:   VEC_ENTRY _c_int04    ;INT04 Externa 
_vector5:   VEC_ENTRY _vec_dummy 
_vector6:   VEC_ENTRY _vec_dummy 
_vector7:   VEC_ENTRY _vec_dummy 
_vector8:   VEC_ENTRY _vec_dummy 
_vector9:   VEC_ENTRY _vec_dummy 
_vector10:  VEC_ENTRY _vec_dummy 
_vector11:  VEC_ENTRY _c_int11      ;ISR address 
_vector12:  VEC_ENTRY _vec_dummy 
_vector13:  VEC_ENTRY _vec_dummy 
_vector14:  VEC_ENTRY _vec_dummy 
_vector15:  VEC_ENTRY _vec_dummy

Biquad bandpass filter bank

● February 2, 2011 ● Coded in Mixed C and ASM for the SHARC

// Efficient implementation of biquad bandpass filter bank.
// This is the core computation for musical physical models ("modal models").
//
// This code takes 1.5 instructions to compute each floating-point bandpass filter;
// for a SHARC running at 3ns instruction cycle this is 4.5ns / biquad.
//
// Written by Lippold Haken of Haken Audio, 2010, using ADSP-21364 and VisualDSP++ C.
//
//  dmXmod - pointer to DM data memory
//  pmXmod - pointer to PM data memory
//  bfbXComputePairs - optimized computation of a bank of bandpass filters
//  bfbXCoef - coefficient computation for the bandpass filters

// CFDSP hardware-specific definitions.
#include "cfdsp21364.h"

// === Start of definitions for .h file
// The following structure definitions are shared by this code and by the calling code,
// so they should be moved to a common .h file.

// Data in dm.
typedef struct 
{
    #define bfb4FilterSections  64      // number of modes (number of bandpass filter sections)
    float bfb_Y[2][bfb4FilterSections];         // save state: y[n-1] and y[n-2]
} DmX_bfb;

// Data in pm.
typedef struct 
{
    // The bfb_C[] contains 3 coefficients for each biquad, with a pair of biquad's coefficients 
    // interleaved so that biquads can be processed in pairs by the optimized SIMD assembly loop.
    float bfb_C[ 3 * bfb4FilterSections ];      // bfb filter coefficients, see comment above
} PmX_bfb;

// Pointers to data, and a samples counter.
DmX_bfb *dmXmod;
PmX_bfb *pmXmod;
int sampInFrame;                    // sample counter

#define sr 48000                    // sample rate in Hz
#define Pi    (3.14159265)
#define ABS(x) ( ((x)>(0)) ? (x) : -(x) )
void sinCos( float fRadians, float *fSin, float *fCos ); // see Lippold Haken's sinCos code snippet 

// === End of definitions for .h file

// Parameter arrays for Biquad Filter Bank (modal physical modelling) synthesis.
//
// If P is the biquad pair number, and R=0/R=1 distinguishes the two biquads within the pair:
//  x[n]   coefficient for biquad k is at pmXmod->bfb_C[ 6*P + R + 0 ]; this is also -x[n-2] coefficient.
//  y[n-1] coefficient for biquad k is at pmXmod->bfb_C[ 6*P + R + 2 ]; this is scaled by 0.5 to avoid overflow
//  y[n-2] coefficient for biquad k is at pmXmod->bfb_C[ 6*P + R + 4 ]
// The first index to dmXmod->bfb_Y[] is based on lsbs of sample counter.
// The second index to dmXmod->bfb_Y[] is 2P+R.
// Note: R=0 always for even voices, R=1 always for odd voices.
#define bfb_Ynm1(P) &(dmXmod->bfb_Y[sampInFrame & 1][2*(P)])        // start in y[n-1][] array
#define bfb_Ynm2(P) &(dmXmod->bfb_Y[1-(sampInFrame & 1)][2*(P)])    // start in y[n-2][] array

void bfbXComputePairs(int biquadP, int pairs, 
            float inDiff0,          // x[n]-x[n-2] for even biquads
            float inDiff1,          // x[n]-x[n-2] for odd biquads
            float *fSum0,           // output for even biquad sum
            float *fSum1)           // output for odd biquad sum
// Sum the output of a sequence of at least 2 biquad pairs.
// For each of the sequence of pairs, the first (even) biquad in each pair has a common input "input0",
// and the second (odd) biquad of each pair has a (possibly different) common input "input1". 
{
    // Assmebly-coded bfb filter loop using floating-point math.
    // The core loop takes 3 instructions for 2 biquads, or 4.5ns per biquad.
    // We use the even-numbered biquads for one voice, the odd-numbered biquads for a second voice.
    asm(
    "#include <def21364.h>                          \n"
    "BIT SET MODE1 SIMD;                            \n"

    // For each biquad k:  (k=6*biquadP for R=0)
    //   y[k](n) = C0[k] * (x[k](n) - x[k](n-2))
    //           + C1[k] * 2 * y[k](n-1) 
    //           + C2[k] * y[k](n-2)
    //
    // Register usage for biquad k (R==0):
    //   F0       = C0[k], C1[k], and C2[k]
    //   F4       = input0-old_input0 aka x[k](n)-x[k](n-2)
    //   F6       = y[k](n-1), y[k](n-2)
    //   F8       = scratch for computing y[k](n)
    //   F10      = sum of even biquad's y[](n) in each biquad pair
    //   F12      = scratch for computing y[k](n)
    //   F13      = final value of y[k-2](n)
    //   DM(I4)   = C0[k], C1[k], and C2[k]
    //   PM(I10)  = read y[k](n-2) write to y[k-2](n)
    //   PM(I12)  = read y[k](n-1)
    // Register usage for biquad k+1 (R==1):
    //   SF0      = C0[k+1], C1[k+1], and C2[k+1]
    //   SF4      = input1-old_input1 aka x[k+1](n)-x[k+1](n-2)
    //   SF6      = y[k+1](n-1), y[k+1](n-2)
    //   SF8      = scratch for computing y[k+1](n)
    //   SF10     = sum of odd biquad's y[](n) in each biquad pair
    //   SF12     = scratch for computing y[k+1](n)
    //   SF13     = final value of y[k-1](n)
    //   DM(I4+1) = C0[k+1], C1[k+1], and C2[k+1]
    //   PM(I10+1)= read y[k+1](n-2) write y[k-1](n)
    //   PM(I12+1)= read y[k+1](n-1)
    // Register usage for even and odd biquads:
    //   M4=M12=2
    //   M10=+4
    //   M11=-2
    //
    // The three instruction loop below computes two biquads simultaneously, for k and k+1,
    // in each SIMD loop execution.  The comments are written just for biquad k;
    // the biquad k+1 code and comments are implied by SIMD mode.

    // Set pointer increments.
    "      M4=2;M10=4;M11=-2;M12=2;                 \n"

    // Loop priming, with k=0 (and SIMD k=1).

        // Initialize F10=SF10=0                    fetch C0[k]         
    "      R10=R10-R10,                             F0=DM(I4,M4);                           \n"

        // C0[k]*(x(n)-x(n-2))                      fetch C1[k]         fetch y[k](n-1)
    "      F8=F0*F4,                                F0=DM(I4,M4),       F6=PM(I12,M12);     \n"
        // C1[k] * y[k](n-1)                        fetch C2[k]         fetch y[k](n-2)
    "      F12=F0*F6,                               F0=DM(I4,M4),       F6=PM(I10,M12);     \n"
        // C2[k] * y[k](n-2)    biquad k sum        fetch C0[k+2]
    "      F12=F0*F6,           F8=F8+F12,          F0=DM(I4,M4);                           \n"

    // Main loop, for k = 2..K-2 by 2 (and SIMD for k = 1..K-1 by 2).
    "DO (PC, endx) UNTIL LCE;   \n"  

        // C0[k]*(x(n)-x(n-2))  F13 is y[k-2](n)    fetch C1[k]         fetch y[k](n-1)
    "      F8=F0*F4,            F13=F8+F12,         F0=DM(I4,M4),       F6=PM(I12,M12);     \n"
        // C1[k] * y[k](n-1)    output sum          fetch C2[k]         fetch y[k](n-2)
    "      F12=F0*F6,           F10=F10+F13,        F0=DM(I4,M4),       F6=PM(I10,M11);     \n"
        // C2[k] * y[k](n-2)    biquad k sum        fetch C0[k+2]       save y[k-2](n)
    "endx: F12=F0*F6,           F8=F8+F12,          F0=DM(I4,M4),       PM(I10,M10)=F13;    \n"

        //                      F13 is y[k-2](n)                                     
    "                           F13=F8+F12,                             modify(I10,M11);    \n"
        //                      output sum                              save y[k-2](n)
    "                           F10=F10+F13,                            PM(I10,M10)=F13;    \n"

    "BIT CLR MODE1 SIMD;                            \n"     // disable SIMD
    "      NOP;                                     \n"     // wait for SIMD disable

    // Output register list, each element: "=rN" (varName)
            :   "=R10" (*fSum0),                    // F10 has floating sum of even & odd biquads
                "=S10" (*fSum1)                     // SF10 has floating sum of even & odd biquads
    // Input register list, each element: "rN" (varName)
            :   "lcntr" (pairs-1),
                "R4" (inDiff0),                     // x[n]-x[n-2] for even biquads
                "S4" (inDiff1),                     // x[n]-x[n-2] for odd biquads
                "I4" (&pmXmod->bfb_C[6 * biquadP]), // coefficients for all biquads
                "I10" (bfb_Ynm2(biquadP)),          // y[n-2] for all biquads, updated to y[n]
                "I12" (bfb_Ynm1(biquadP))           // y[n-1] for all biquads
    // Clobbered register list:  
    // We must avoid "do not use" regs, Compiler and Library Manual p.1-245.
    // We try to use mostly scratch regs, Compiler and Library Manual p.1-248.
            :   "r0", "r4", "r6", "r8", "r10", "r12", "r13", 
                "i4", "i10", "i12", 
                "m4", "m10", "m11", "m12" );
}

void bfbXCoef(int biquadP, int biquadR, float freq, float amp, float bw)
// Compute coefficients for a biquad filter section of the floating-point Biquad Filter Bank routine.
{
    int startingIndex = 6 * biquadP + biquadR;          // 6 * biquad pair number (+1 if second in pair)
    float *pC = &pmXmod->bfb_C[ startingIndex ];        // index biquad filter coefficients array
    float *y0 = &dmXmod->bfb_Y[0][2*biquadP+biquadR];   // filter memory for biquad
    float *y1 = &dmXmod->bfb_Y[1][2*biquadP+biquadR];   // filter memory for biquad

    // Is center frequency is below aliasing-cutoff limit?
    if (freq != 0.0 && freq < 0.5 * sr)
    {
        // Compute sin and cos of the mode frequency.
        float phase = freq * 2. * Pi / sr;
        float fSin, fCos;
        sinCos(phase, &fSin, &fCos); // efficient sine and cosine; see Lippold Haken's DSPrelated code snippet.

        // These formulas are adapted from Robert Bristow-Johnson BLT biquad web posting.
        // I use the BPF with "constant skirt gain, peak gain = Q", with amplitude scaling of input coefficients.
        // Compute intermediate parameters, alpha and beta.
        //      alpha = sin(w0)/(2*Q)
        //      beta  = 1/(1+alpha)
        //      y[n] = (beta * sin(w0)/2   )  *  (x[n] - x[n-2])
        //           + (beta * 2*cos(w0)   )  *  y[n-1]
        //           + (beta * (alpha - 1) )  *  y[n-2] 
        // In addition, the input x-coefficients are scaled by the amplitude of the biquad --
        // I use **twice** the amplitude value.
        float alpha = fSin * bw / 2.0;      // this is sin/(2*q)
        float beta = 1.0 / (1.0 + alpha); 
        pC[0] = amp * beta * fSin;          // x[n] coefficient, and -x[n-2] coefficient
        pC[2] = beta * 2. * fCos;           // y[n-1] coefficient
        pC[4] = beta * (alpha - 1.0);       // y[n-2] coefficient
    }
    else
    {
        // Frequency of bfb filter is beyond aliasing frequency,
        // or if we are about to start a new note.
        // Zero bfb filter coefficients and zero out the filter memory.
        pC[0] = pC[2] = pC[4] = 0;
        *y0 = *y1 = 0;
    }
}

2D IFFT

Senthilkumar ● February 2, 2011 ● Coded in Scilab

function [a] =ifft2d(a2) 
//a2 = 2D-DFT of any real or complex 2D matrix
//a = 2D-IDFT of a2  
m=size(a2,1)
n=size(a2,2)
//Inverse Fourier transform along the rows
for i=1:n
a1(:,i)=exp(2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a2(:,i) 
end
//Inverse fourier transform along the columns
for j=1:m
atemp=exp(2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).' 
a(j,:)=atemp.'
end
a = a/(m*n)
a = real(a)
endfunction

2D FFT

Senthilkumar ● February 2, 2011 ● Coded in Scilab

function [a2] = fft2d(a) 
//a = any real or complex 2D matrix
//a2 = 2D-DFT of 2D matrix  'a'
m=size(a,1)
n=size(a,2)
// fourier transform along the rows
for i=1:n
a1(:,i)=exp(-2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a(:,i) 
end
// fourier transform along the columns
for j=1:m
a2temp=exp(-2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).' 
a2(j,:)=a2temp.'
end
for i = 1:m
    for j = 1:n
        if((abs(real(a2(i,j)))<0.0001)&(abs(imag(a2(i,j)))<0.0001))
            a2(i,j)=0;
        elseif(abs(real(a2(i,j)))<0.0001)
            a2(i,j)= 0+%i*imag(a2(i,j));
        elseif(abs(imag(a2(i,j)))<0.0001)
            a2(i,j)= real(a2(i,j))+0;
        end
    end
end

Circular Convolution Matrix

● February 1, 2011 ● Coded in Matlab

function HC = convmtx_circ(h, N);
% CONVMTX_CIRC Circular Convolution Matrix.
%    CONVMTX_CIRC(h,N) returns the circular convolution matrix for vector h.
%    h must be a column vector. If X is a column vector of length N,
%    then CONVMTX_CIRC(h,N)*X is the same as the circular convolution h * X.
%
% by Danilo Zanatta - 01.02.2011

h = h(:);
Nh = length(h);
Nc = N + Nh - 1;

if N < Nh,
    error('The circular convolution matrix size must be greater or equal the channel length');
end

H = convmtx(h,N);

HC = H(1:N, 1:N);
HC(1:(Nh-1), 1:N) = HC(1:(Nh-1), 1:N) + H((N+1):Nc, 1:N);

return

This function helps in converting linear phase fir filter transfer function to min phase fir filter transfer function

● January 31, 2011●1 comment ● Coded in Matlab

function [hn_min] = lp_fir_2_mp_fir(hn_lin)

% This function helps in converting the linear phase 
% FIR filter to minimum phase FIR filter. It essentially 
% brings all zeros which are outside the unit circle to 
% inside of unit circle.

r 	= roots(hn_lin);
k	= abs(r) > 1.01;
r1  	= r(k);		% roots outside of unit circle.
r2	= r(~k);% roots on or within the unit circle
s	= [1./r1; r2];
hn_min	= poly(s);
hn_min	= hn_min*sqrt(sum(hn_lin.^2))/sqrt(sum(hn_min.^2));

GPIO.c - Using the GPIO as output

David Valencia ● January 31, 2011●4 comments ● Coded in C for the TI C67x

// GPIO.c	UPIITA-IPN 
// JOSE DAVID VALENCIA PESQUEIRA
//
// THIS PROGRAM ALLOWS THE DSP TO SEND A BIT THROUGH A GPIO PIN (PIN 7 IN THIS EXAMPLE)
// TO TURN ON A LED ON

#include <stdio.h>
#include <stdlib.h>
#include <csl_gpio.h>
#include <csl_gpiohal.h>
#include <csl_irq.h>

// GLOBAL VARIABLES
volatile int flag = 1;
volatile int pulse_flag = 1;

// CODE TO DEFINE GPIO HANDLE
GPIO_Handle gpio_handle; 

// GPIO REGISTER CONFIGURATION
GPIO_Config gpio_config = {         
    0x00000000, // gpgc = Interruption passthrough mode
    0x0000FFFF, // gpen = All GPIO 0-15 pins enabled
    0x0000FFFF, // gdir = All GPIO pins as outputs
    0x00000000, // gpval = Saves the logical states of pins
    0x00000000, // gphm All interrupts disabled for IO pins
    0x00000000, // gplm All interrupts for CPU EDMA disabled
    0x00000000  // gppol -- default state */
};

// Function prototypes
void send_edge();
void delay();

// Function definitions
void delay()
{
	// Delay function
	int count, count2;
	for(count = 0; count < 200; count++)
	{
		for(count2 = 0; count2 < 50000; count2++);
	}
}

void send_edge()
{	
    DSK6713_init();  
    // Open and configure the GPIO
    gpio_handle = GPIO_open( GPIO_DEV0, GPIO_OPEN_RESET );
    GPIO_config(gpio_handle,&gpio_config);
    // Send values through the pins
	
    GPIO_pinWrite(gpio_handle,GPIO_PIN7,0);
    delay();
    GPIO_pinWrite(gpio_handle,GPIO_PIN7,1);
    delay();
    GPIO_pinWrite(gpio_handle,GPIO_PIN7,0);
}

main()
{     
	send_edge(); // Configures the pins and sends a rising edge
	printf(“ ¡Felicidades! ¡Has logrado encender el led! ”);
	
}
// End of program>>

Inverse Sqrt and Fourth Root approximations

● January 31, 2011 ● Coded in C for the SHARC

// Fast InvSqrt approxumation, with an error of less than 4%. 
// This approximation is attributed to Greg Walsh.
inline void InvSqrt(float *x) {	*(int *)x = 0x5f3759df - (*(int *)x >> 1); } 
inline float SqrtSqrt(float x) { InvSqrt(&x); InvSqrt(&x); return x; }

Fourier Trigonometric Series Approximation

David Valencia ● January 31, 2011 ● Coded in Matlab

%Trigonometric Fourier Series Approximation
%José David Valencia Pesqueira - UPIITA-IPN
%As posted for Dsprelated.com
 
close all; clear all; clc;
syms t k;

fprintf('\n ### FOURIER TRIGONOMETRICAL SERIES###');

fprintf('\nSpecify the range (t1,t2) for the approximation:\n');
t1=input('t1= ');
t2=input('t2= ');

% Getting the period

fprintf('\nIs there a repeating patterin in this range? (Y/N): ');
resp1=input('','s');

switch resp1
    case 'Y'
        fprintf('\nWhat is the period for the function?\n');
        T=input('T= ');
        omega0=(2*pi)/T;
    case 'N'
        T=t2-t1;
        omega0=(2*pi)/T;
    otherwise
        fprintf('\nInvalid Option');
        fprintf('\nEnd of program>>>');
        return
end

fprintf('\nThe approximation will be generated in a trigonometrical series: \n\n');
fprintf('f(t)=a0+a1*cos(omega0*t)+a2*cos(2*omega0*t)+...+b1*sen(omega0*t)+b2*sen(2*omega0*t)');
 
fprintf('\n\nUntil what coefficient ak,bk should I compute? (positive integer) n= ');
n=input('');
if n<=0
    fprintf('Invalid n');
    fprintf('\nProgram Stop >>>');
    return
end

fprintf('\nIs the function defined by parts? (Y/N): ');
resp1=input('','s');

switch resp1
    case 'N'
        fprintf('\nWrite the function in terms of t\nf(t)= ');
        f=input('');
        
        faprox=0;

        a0=(1/T)*int(f,t,t1,(t1+T));
 
        faprox=a0;
 
        for k=1:n
            a(k)=(2/T)*int((f*cos(k*omega0*t)),t,t1,(t1+T));
            faprox=faprox+a(k)*cos(k*omega0*t);
        end
 
        for k=1:n
            b(k)=(2/T)*int((f*sin(k*omega0*t)),t,t1,(t1+T));
            faprox=faprox+b(k)*sin(k*omega0*t);
        end

    % If the function is defined in many parts
    case 'Y'
        fprintf('\nHow many ranges do you want to define (positive integer): ');
        numinterv=input('');
        numinterv=round(numinterv);
        for k=1:numinterv
            fprintf('\n## Range #%d definition',k);
            fprintf('\n Write the range in the form of t(%d)<t<t(%d) : ',k-1,k);
            fprintf('\nt(%d)= ',k-1);
            a(k,1)=input(''); %Initial limits vector
            fprintf('\nt(%d)= ',k);
            a(k,2)=input(''); %Final limit vector
            fprintf('\nThe range # %d has a start in %d and ends in %d',k,a(k),a(k+1));
            fprintf('\nPlease define f(t) for the %d th interval:\nf(t)= ',k);
            ft(k)=input('');
            fprintf('\nEnd of definition of the function f%d',k);
        end
        
        faprox=0;
        a0=0
        for j=1:numinterv
            a0=a0+(1/T)*int(ft(j),t,a(j,1),a(j,2));
        end
 
        faprox=a0;
 
        % ak coefficient
        ak=0;
        for j=1:numinterv
            acumulador=(2/T)*int((f(j)*cos(k*omega0*t)),t,a(j,1),a(j,2));
            ak=ak+acumulador;
        end

        return
        
        fprintf('\n## End of Debug Execution');
        return    
    otherwise
        fprintf('\nEND OF PROGRAM>>');
        return
end

%% Salida de datos
fprintf('\nThe resulting coeficcients are: ');
a0
a
b
 
fprintf('\n Do you wish to plot the graph of f(t) against the approximate series? (Y/N): ');
resp1=input('','s');
 
if resp1=='Y'
    ezplot(f,[t1,t2]);
    hold on
    ezplot(faprox,[t1,t2]);
    grid;
end
 
fprintf('\nEND OF PROGRAM');

Fast power-of-10 approximation, for 32-bit floats