## Delay estimation revisited

June 9, 20125 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

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

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';

% *************************************************
% *************************************************
switch smode
case 'evalConventional'
% conventionally designed RRC filter

case 'evalIdeal'
h0 = 1;

case 'evalEqualized'
% alternative RRC design that equalizes the known frequency
% response of the resampler

otherwise assert(false);
end

% *************************************************
% --- 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, , 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, , 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, , 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, , 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, , 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

December 18, 20114 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);

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

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

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

December 10, 20111 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

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

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 %.
% *************************************************************
% 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
% *************************************************************
% ---------------------------
% The baseband-equivalent signal represents the modulation on
% 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

sig (1) = A;

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

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

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

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

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

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

% exclude carrier:
sig(1) = 0;

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

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

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

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

## GSM (GMSK) power spectrum equation

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.
% 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];

xlabel('f/MHz');
ylabel('PSD/dB');
title('GSM spectrum');``````

## Yet Another FIR design algorithm

August 21, 20113 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'
% ********************************************
% ********************************************
% 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;

% sanity check... (orthogonality of wanted and unwanted component)
assert(sum(mask_inband) > 0, 'passband is empty');
'passband and stopband overlap');

% ********************************************
% start with flat passband signals at input and output of filter
% those will become the input to the LMS solver.
% ********************************************
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
% => 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

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

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
% ********************************************
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 = doPlot_dB(d, EVM_dB, 'r', 'error');

handle = plot(d.fb(inbandBins), d.errorSpecBinVal_inband_dB(inbandBins), 'r+');
set(handle, 'markersize', 2);

d=doPlotEnd(d);
ylim([-100, 10]);

if false
% ********************************************
% plot constraint signal and weights
% ********************************************
figure(31); grid on; hold on;
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);
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);
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

% region where second argument is defined

% region where result will be defined
switch(p.logicalOperator)
case 'or'
case 'and'
otherwise
assert(false);
end

% 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);
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

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