DSPRelated.com
Code
   Code Snippets    Submitted by mnentwig

Delay estimation revisited

Markus Nentwig ● June 9, 2012●5 comments ● Coded in Matlab 
% ****************************************************************
% find sub-sample delay and scaling factor between two cyclic signals 
% to maximize crosscorrelation
% Markus Nentwig, 120609_v1.1
% ****************************************************************
function iterDelayEstDemo();
    close all;
    n = 1024;    
    % n = 1024 * 256; disp('*** test: long signal enabled ***');

    % ****************************************************************
    % random signal
    % ****************************************************************
    fd = randn(1, n) + 1i * randn(1, n);
    
    % ****************************************************************
    % lowpass filter
    % ****************************************************************
    f = binFreq(n);
    fd(abs(f) > 0.045) = 0;
    s1 = real(ifft(fd)) * sqrt(n);
    
    % ****************************************************************
    % create delayed 2nd signal
    % ****************************************************************
    dTest_samples = 12.3456;
    cTest = 1.23456;
    % cTest = cTest + 1i; disp('*** test: complex coeff enabled ***'); 
    % cTest = -cTest; disp('*** test: negative coeff enabled ***'); 

    s2 = cTest * cs_delay(s1, 1, dTest_samples);
    % s2 = s2 + 0.5*randn(size(s2)); disp('*** test: noise enabled ***');
    
    % ****************************************************************
    % estimate delay
    % ****************************************************************
    [delay_samples, coeff] = iterDelayEst(s1, s2);

    % ****************************************************************
    % correct it
    % ****************************************************************
    s2a = cs_delay(s1, 1, delay_samples);
    s2b = s2a * coeff;

    figure(); hold on;
    h = plot(real(s1), 'k'); set(h, 'lineWidth', 3);
    h = plot(real(s2), 'b'); set(h, 'lineWidth', 3);
    h = plot(real(s2a), 'r'); set(h, 'lineWidth', 1);
    h = plot(real(s2b), 'm'); set(h, 'lineWidth', 1);
    xlim([1, numel(s1)]);
    xlabel('samples');
    legend('s1', 's2', 's2 un-delayed', 's2 un-delayed and scaled');
    title('test signals');
    
    format long;
    disp('nominal delay of s2 relative to s1')';
    dTest_samples
    disp('iterDelayEst() returned:');
    delay_samples
    disp('original scaling factor:');
    cTest
    disp('estimated scaling factor:');
    coeff
end

% ****************************************************************
% estimates delay and scaling factor 
% ****************************************************************
function [delay_samples, coeff] = iterDelayEst(s1, s2)
    
    s1 = s1(:) .'; % force row vectors
    s2 = s2(:) .';
    rflag = isreal(s1) && isreal(s2);
    n = numel(s1);
    halfN = floor(n/2);
    assert(numel(s2) == n, 'signals must have same length');

    % ****************************************************************
    % constants
    % ****************************************************************    
    % exit if uncertainty below threshold
    thr_samples = 1e-7;

    % exit after fixed number of iterations
    nIter = 25;

    % frequency domain representation of signals
    fd1 = fft(s1);
    fd2 = fft(s2);    

    % first round: No delay was applied
    tau = [];
    fd2Tau = fd2; % delayed s2 in freq. domain
    
    % frequency corresponding to each FFT bin -0.5..0.5
    f = binFreq(n);

    % uncertainty plot data
    e = [];

    % normalization factor
    nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n; % normalizes to 1
    
    % search window: 
    % known maximum and two surrounding points
    x1 = -1;
    x2 = -1;
    x3 = -1;
    y1 = -1;
    y2 = -1;
    y3 = -1;
    
    % ****************************************************************
    % iteration loop
    % ****************************************************************
    for count = 1:nIter
    
        % ****************************************************************
        % crosscorrelation with time-shifted signal
        % ****************************************************************
        xcorr = abs(ifft(fd2Tau .* conj(fd1)))/ nf;

        % ****************************************************************
        % detect peak
        % ****************************************************************
        if isempty(tau)

            % ****************************************************************
            % startup
            % initialize with three adjacent bins around peak
            % ****************************************************************
            ix = find(xcorr == max(xcorr));
            ix = ix(1); % use any, if multiple bitwise identical peaks
            
            % indices of three bins around peak
            ixLow = mod(ix-1-1, n) + 1; % one below
            ixMid = ix;
            ixHigh = mod(ix-1+1, n) + 1; % one above

            % delay corresponding to the three bins
            tauLow = mod(ixLow -1 + halfN, n) - halfN;
            tauMid = mod(ixMid -1 + halfN, n) - halfN;         
            tauHigh = mod(ixHigh -1 + halfN, n) - halfN; 

            % crosscorrelation value for the three bins
            xcLow = xcorr(ixLow);
            xcMid = xcorr(ixMid);
            xcHigh = xcorr(ixHigh);
            
            x1 = tauLow;
            x2 = tauMid;
            x3 = tauHigh;
            y1 = xcLow;
            y2 = xcMid; 
            y3 = xcHigh;
        else
            % ****************************************************************
            % only main peak at first bin is of interest
            % ****************************************************************
            tauMid = tau;
            xcMid = xcorr(1);

            if xcMid > y2
                % ****************************************************************
                % improve midpoint
                % ****************************************************************
                if tauMid > x2
                    % midpoint becomes lower point
                    x1 = x2;
                    y1 = y2;
                else
                    % midpoint becomes upper point
                    x3 = x2;
                    y3 = y2;
                end
                x2 = tauMid;
                y2 = xcMid;
            
            elseif tauMid < x2
                % ****************************************************************
                % improve low point
                % ****************************************************************
                assert(tauMid >= x1); % bitwise identical is OK
                assert(tauMid > x1 || xcMid > y1); % expect improvement
                x1 = tauMid;
                y1 = xcMid;
            elseif tauMid > x2 
                % ****************************************************************
                % improve high point
                % ****************************************************************
                assert(tauMid <= x3); % bitwise identical is OK                
                assert((tauMid < x3) || (xcMid > y3)); % expect improvement
                x3 = tauMid;
                y3 = xcMid;
            else
                assert(false, '?? evaluated for existing tau ??');
            end
        end

        % ****************************************************************
        % calculate uncertainty (window width)
        % ****************************************************************
        eIter = abs(x3 - x1);
        e = [e, eIter];
        if eIter < thr_samples
            % disp('threshold reached, exiting');
            break;
        end

        if y1 == y2 || y2 == y3
            % reached limit of numerical accuracy on one side
            usePoly = 0;
        else
            % ****************************************************************
            % fit 2nd order polynomial and find maximum
            % ****************************************************************
            num = (x2^2-x1^2)*y3+(x1^2-x3^2)*y2+(x3^2-x2^2)*y1;
            denom = (2*x2-2*x1)*y3+(2*x1-2*x3)*y2+(2*x3-2*x2)*y1;
            if denom ~= 0
                tau = num / denom;
                % is the point within [x1, x3]? 
                usePoly = ((tau > x1) && (tau < x3));
            else
                usePoly = 0;
            end            
        end
        if ~usePoly
            % revert to linear interpolation on the side with the
            % less-accurate outer sample 
            % Note: There is no guarantee that the side with the more accurate
            % outer sample is the right one, as the samples aren't 
            % placed on a regular grid!
            % Therefore, iterate to improve the "worse" side, which will
            % eventually become the "better side", and iteration converges.
            
            tauLow = (x1 + x2) / 2;
            tauHigh = (x2 + x3) / 2;
            if y1 < y3
                o = [tauLow, tauHigh];
            else
                o = [tauHigh, tauLow];                
            end
            % don't choose point that is identical to one that is already known
            tau = o(1);
            if tau == x1 || tau == x2 || tau == x3
                tau = o(2);
                if tau == x1 || tau == x2 || tau == x3
                    break;
                end
            end
        end

        % ****************************************************************
        % advance 2nd signal according to location of maximum
        % phase shift in frequency domain - delay in time domain
        % ****************************************************************
        fd2Tau = fd2 .* exp(2i * pi * f * tau);
    end % for

    % ****************************************************************
    % plot the uncertainty (window size) over the number of iterations
    % ****************************************************************
    if true
        figure(); semilogy(e, '+-'); grid on;
        xlabel('iteration');
        title('uncertainty in delay');
    end

    % ****************************************************************
    % the delay estimate is the final location of the delay that 
    % maximized crosscorrelation (center of window).
    % ****************************************************************
    delay_samples = x2;

    % ****************************************************************
    % Coefficient: Turn signal 1 into signal 2
    % ****************************************************************
    coeff = fd2Tau * fd1' ./ (fd1 * fd1');

    % ****************************************************************
    % chop roundoff error, if input signals are known to be 
    % real-valued.
    % ****************************************************************
    if rflag
        coeff = real(coeff);
    end
end

% ****************************************************************
% frequency corresponding to FFT bin
% ****************************************************************
function f = binFreq(n)
    f = (mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
end

% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
    rflag = isreal(waveform);
    
    n = numel(waveform);
    cycLen_s = n / rate_Hz;
    nCyc = delay_s / cycLen_s();
    
    f = 0:(n - 1);
    f = f + floor(n / 2);
    f = mod(f, n);
    f = f - floor(n / 2);
    
    phase = -2 * pi * f * nCyc;
    rot = exp(1i*phase);
    
    waveform = ifft(fft(waveform) .* rot);
    if rflag
        waveform = real(waveform);
    end
end

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

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

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

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

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

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

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

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

Alias/error simulation of interpolating RRC filter

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    if forceReal
        shiftedRef = real(shiftedRef);
    end
end

Baseband-equivalent phase noise model

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

Digital model for analog filters

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Spectrum analyzer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

integrate RMS phase error from spectrum

Markus Nentwig ● October 21, 2011●1 comment ● Coded in Matlab 
% ****************************************************************
% RMS phase error from phase noise spectrum
% Markus Nentwig, 2011
%
% - integrates RMS phase error from a phase noise power spectrum
% - generates an example signal and determines the RMS-average
%   phase error (alternative method)
% ****************************************************************
function pn_snippet()
    close all;
    
    % ****************************************************************
    % Phase noise spectrum model
    % --------------------------
    % Notes: 
    % - Generally, a phase noise spectrum tends to follow
    %   |H(f)| = c0 + c1 f^n1 + c2 f^n2 + c3 f^n3 + ...
    %   Therefore, linear interpolation should be performed in dB
    %   on a LOGARITHMIC frequency axis.
    % - Single-/double sideband definition: 
    %   The phase noise model is defined for -inf <= f <= inf
    %   in other words, it contributes the given noise density at
    %   both positive AND negative frequencies. 
    %   Assuming a symmetrical PN spectrum, the model is evaluated 
    %   for |f| => no need to explicitly write out the negative frequency
    %   side.
    % ****************************************************************
    
    % PN model
    % first column: 
    %   frequency offset
    % second column: 
    %   spectral phase noise density, dB relative to carrier in a 1 Hz 
    %   observation bandwidth

    d = [0 -80
         10e3 -80
         1e6 -140
         9e9 -140];    
    f_Hz = d(:, 1) .';
    g_dBc1Hz = d(:, 2) .' -3;
    
    % get RMS phase error by integrating the power spectrum
    % (alternative 1)
    e1_degRMS = integratePN_degRMS(f_Hz, g_dBc1Hz)

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

% ****************************************************************
% Integrate the RMS phase error from the power spectrum
% ****************************************************************
function r = integratePN_degRMS(f1_Hz, g1_dBc1Hz)
    1;
    % integration step size
    deltaF_Hz = 100;
    
    % list of integration frequencies
    f_Hz = deltaF_Hz:deltaF_Hz:5e6;
    
    % interpolate spectrum on logarithmic frequency, dB scale
    % unit is dBc in 1 Hz, relative to a unity carrier
    fr_dB = interp1(log(f1_Hz+eps), g1_dBc1Hz, log(f_Hz+eps), 'linear');
    
    % convert to power in 1 Hz, relative to a unity carrier
    fr_pwr = 10 .^ (fr_dB/10);

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

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

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

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

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

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

    % indices of positive-frequency FFT bins. 
    % Does not include DC
    ixPos = 2:(n/2);
    
    % indices of negative-frequency FFT bins. 
    % Does not include the bin on the Nyquist limit
    assert(mod(n, 2) == 0, 'n must be even');
    ixNeg = n - ixPos + 2;
    
    % Generate signal in the frequency domain
    sig = zeros(1, n);
    % positive frequencies
    sig(ixPos) = randn(size(ixPos)) + 1i * randn(size(ixPos));
    % for purely real-valued signal: conj()
    % for purely imaginary-valued signal such as phase noise: -conj()
    % Note: 
    %   Small-signals are assumed. For higher "modulation indices", 
    %   phase noise would become a circle in the complex plane
    sig(ixNeg) = -conj(sig(ixPos));
    
    % normalize energy to unity amplitude A per bin
    sig = sig * A / RMS(sig);

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

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

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

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

    % plot on linear f axis
    figure(); hold on;
    plot(fftshift(fb_Hz), fftshift(mag2dB(H)), 'k');
    xlabel('f/Hz'); ylabel('dBc in 1 Hz');
    title('phase noise spectrum (linear frequency axis)');
    
    % plot on logarithmic f axis
    figure(); hold on;
    ix = find(fb_Hz > 0);
    semilogx(fb_Hz(ix), mag2dB(H(ix)), 'k');
    xlabel('f/Hz'); ylabel('dBc in 1 Hz');
    title('phase noise spectrum (logarithmic frequency axis)');
    
    % apply frequency response to signal
    sig = sig .* H;
    
    % add carrier
    sig (1) = A;

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

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

    % RMS average
    % note: exp(1i*x) ~ x
    %       as with magnitude, power/energy is proportional to phase error squared
    r = RMS(ph_deg);
    
    % knowing that the signal does not have an amplitude component, we can alternatively
    % determine the RMS phase error from the power of the phase noise component

    % exclude carrier:
    sig(1) = 0; 

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

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

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

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

GSM (GMSK) power spectrum equation

Markus Nentwig ● September 18, 2011 ● Coded in Matlab 
% ************************************************************
% Spectrum model for GSM signal 
% Markus Nentwig, 2011
% based on 3GPP TS 45.004 section 2 "Modulation format for GMSK", 
% assuming no additional filtering and an infinite-length
% symbol stream (no slot structure)
% ************************************************************

% ************************************************************
% Parameters
% The "baseline" serves as a gentle reminder that the model
% is only valid near the center frequency.
% For large frequency offsets, one would need more information on
% the particular transmitter used.
% If not available, assume that spectral emission mask requirements
% are met.
% ************************************************************
BW = 285e3;
BW2 = 186e3;
baseline_dB = -76;
% baseline_dB = -999 % disable the constant term

% ************************************************************
% empirical GSM (GMSK narrow-bandwidth pulse) model equation
% ************************************************************
f = (-500e3:1e3:500e3)+1e-3;
gaussPart = exp(-(2*f/BW) .^2);
sincPart = sin(pi*f/BW2) ./ (pi*f/BW2);
flatPart = 10^(baseline_dB/20);
H_dB = 20*log10(abs(gaussPart .* sincPart) + flatPart);

% ************************************************************
% plot the spectrum
% ************************************************************
figure(); grid on; hold on;
h = plot(f/1e6, H_dB, 'b'); set(h, 'linewidth', 2);

% ************************************************************
% plot 'a' GSM spectral emission mask
% note, this is only "an" example
% See 3GPP TS 45.005 section 4.2.1
%   "Spectrum due to the modulation and wide band noise"
% section 4.2.1.1 
%   "General requirements for all types of Base stations and MS"
% note the warning regarding measuring the nominal signal level!
% ************************************************************
maskX_MHz = [0, 100e3, 200e3, 250e3, 400e3, 600e3]/1e6;
maskY_dB = [0.5, 0.5, -30, -33, -60, -60];
h = plot(maskX_MHz, maskY_dB, 'k'); set(h, 'linewidth', 3);

legend('|H(f)|', 'GSM mask example');
xlabel('f/MHz');
ylabel('PSD/dB');
title('GSM spectrum');

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