## Function to perform Adaptive filtering

December 28, 2011 Coded in Scilab
``````function[y,M] = adapt_filt(xlms,B,h,delta,l,l1)
//x = the signal from the speaker directly
//B = the signal thorugh the echo path
//h = impulse response of adaptive filter
//l = length of signal 'x'
//l1 = length of adaptive filter order
for k = 1:150
for n = 1:l
xcap = xlms((n+l1-1):-1:(n+l1-1)-l1+1);
yout(n) = h*xcap;
e(n) = B(n)- yout(n);
xnorm = 0.001+(xcap*xcap');
h = h+((delta*e(n))*(xcap'));
end
eold = 0.0;
for i = 1:l
MSE = eold+(e(i)^2);
eold = MSE;
end
if MSE <= 0.0001 then
break;
end
end
y =  zeros(1,length(e));
M = zeros(1,length(h));
y  = e;
M = h;
endfunction``````

## Full band Echo canceller

December 28, 2011 Coded in Scilab
``````//Caption : FULL Band Echo Cnacellation using NLMS Adaptive Filter

clc;
order = 40;  // Adaptive filter order
x = x';
N = length(x);  //length of speech signal
//Delay introduced in echo path
delay = 100;
//Echo at same speaker
xdelayed = zeros(1,N+delay);
for i = delay+1:N+delay
xdelayed(i) = x(i-delay);
end
//Initialize the adaptive filter coefficients to zero
hcap = zeros(1,order);
//To avoid negative values generated during convolution
for i = 1:order-1
xlms(i) = 0.0;
end
for i = order:N+order-1
xlms(i) = x(i-order+1);
end
Power_X = pow_1(x,N);  //Average power of speech signal
//Calculation of step size of adaptive filter
delta = 1/(10*order*Power_X);
figure(1)
subplot(3,1,1)
plot([1:N],x)
title('Speech Signal Generated by some Speaker A')
sound(x,Fs,16)
subplot(3,1,2)
plot([1:length(xdelayed)],xdelayed)
title('Echo signal of speaker A received by speaker A')
sound(xdelayed,Fs,16)
subplot(3,1,3)
plot([1:length(x_out)],x_out)
title('Echo signal at speaker A after using NLMS Adaptive filter Echo Canceller')
figure(2)
title('Adaptive Filter (Echo Canceller) Impulse response')``````

## Chaning the bit depth in speech samples

December 28, 2011 Coded in Scilab
``````//Caption: Reading a Speech Signal &
//. Write it in another file
//. Changing the bit depth from 16 bits to 8 bits
clear;
clc;
wavwrite(y,Fs,8,"E:\4_8bit.wav");
Ny = length(y); //Number of samples in y (4.wav)
Nx = length(x); //Number of samples in x (4_8bit.wav)
Memory_y = Ny*bits_y;  //memory requirement for 4.wav in bits
Memory_x = Nx*bits_x;  //memory requirement for 4_8bit.wav in bits
disp(Memory_y,'memory requirement for 4.wav in bits =')
disp(Memory_x,'memory requirement for 4_8bit.wav in bits =')
//Result
//memory requirement for 4.wav in bits =
//     133760.
//
//memory requirement for 4_8bit.wav in bits =
//     66880.``````

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

## Python zplane function

December 17, 20114 comments Coded in Python
``````#
# Copyright (c) 2011 Christopher Felton
#
# This program is free software: you can redistribute it and/or modify
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# The following is derived from the slides presented by
# Alexander Kain for CS506/606 "Special Topics: Speech Signal Processing"
# CSLU / OHSU, Spring Term 2011.

import numpy as np
import matplotlib.pyplot as plt
from  matplotlib import patches
from matplotlib.figure import Figure
from matplotlib import rcParams

def zplane(b,a,filename=None):
"""Plot the complex z-plane given a transfer function.
"""

# get a figure/plot
ax = plt.subplot(111)

# create the unit circle
color='black', ls='dashed')

# The coefficients are less than 1, normalize the coeficients
if np.max(b) > 1:
kn = np.max(b)
b = b/float(kn)
else:
kn = 1

if np.max(a) > 1:
kd = np.max(a)
a = a/float(kd)
else:
kd = 1

# Get the poles and zeros
p = np.roots(a)
z = np.roots(b)
k = kn/float(kd)

# Plot the zeros and set marker properties
t1 = plt.plot(z.real, z.imag, 'go', ms=10)
plt.setp( t1, markersize=10.0, markeredgewidth=1.0,
markeredgecolor='k', markerfacecolor='g')

# Plot the poles and set marker properties
t2 = plt.plot(p.real, p.imag, 'rx', ms=10)
plt.setp( t2, markersize=12.0, markeredgewidth=3.0,
markeredgecolor='r', markerfacecolor='r')

ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

# set the ticks
r = 1.5; plt.axis('scaled'); plt.axis([-r, r, -r, r])
ticks = [-1, -.5, .5, 1]; plt.xticks(ticks); plt.yticks(ticks)

if filename is None:
plt.show()
else:
plt.savefig(filename)

return z, p, k``````

## Testing the Flat-Top Windowing Function

December 14, 2011 Coded in Matlab
``````%    Code for testing the 'Wind_Flattop(Spec)' function in
%    reducing 'scalloping loss' errors in time signal amplitude
%    estimation.
%
%    Generates a time-domain sinusoid, computes its FFT, and passes
%    that FFT sequence to the 'Wind_Flattop(Spec)' function.
%
%    The maximum output sample of the 'Wind_Flattop(Spec)' function
%    is used to estimate the peak amplitude of the original
%    sinusoidal time-domain test signal.
%
%    The user controls the values for the test sinusoid's
%    'Test_Freq' and peak amplitude 'Peak_Amp' near lines 23 & 24.
%
%    Richard Lyons [December, 2011]

clear all, clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Define test parameters
Test_Freq = 7.22; % Test tone's freq. Must be less that N/2
Peak_Amp = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N = 64; % Number of time samples
Index = (0:N-1);

X = Peak_Amp*cos(2*pi*(Test_Freq)*Index/N + pi/3);

figure(1), clf
subplot(2,1,1)
plot(X,':ko', 'markersize', 4)
title('Original time signal'), grid on, zoom on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   FFT the input 'X' sequence and call 'Wind_Flattop()' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Spec = fft(X);
[Windowed_Spec] = Wind_Flattop(Spec);

subplot(2,1,2), plot(abs(Spec),'ko', 'markersize', 3)
title('SpecMag of unwindowed time signal'), grid on, zoom on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Display results accuracy (error)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(['Test Freq = ',num2str(Test_Freq),...
',  True Peak Amplitude = ',num2str(Peak_Amp)])

Mag_peak_unwindowed = max(abs(Spec));
Unwindowed_Amp_Estimate = 2*Mag_peak_unwindowed/N;
Unwindowed_Amp_Estimate_Error_in_dB = ...
20*log10(Unwindowed_Amp_Estimate/Peak_Amp);
disp(' ')
disp(['Unwindowed Peak Amplitude Estimate = ',...
num2str(Unwindowed_Amp_Estimate)])
disp(['Unwindowed Estimate Error in dB = ',...
num2str(Unwindowed_Amp_Estimate_Error_in_dB),' dB'])

M_peak_windowed = max(abs(Windowed_Spec));
Windowed_Amp_Estimated = 2*M_peak_windowed/N;
Windowed_Amp_Estimation_Error_in_dB = ...
20*log10(Windowed_Amp_Estimated/Peak_Amp);
disp(' ')
disp(['Windowed Peak Amplitude Estimate = ',...
num2str(Windowed_Amp_Estimated)])
disp(['Windowed Estimate Error in dB = ',...
num2str(Windowed_Amp_Estimation_Error_in_dB),' dB'])``````

## Flat-Top Windowing Function for the Accurate Measurement of a Sinusoid's Peak Amplitude Based on FFT Data

December 14, 2011 Coded in Matlab
``````function [Windowed_Spec] = Wind_Flattop(Spec)

%  Given an input spectral sequence 'Spec', that is the
%  FFT of some time sequence 'x', Wind_Flattop(Spec)
%  returns a spectral sequence that is equivalent
%  to the FFT of a flat-top windowed version of time
%  sequence 'x'.  The peak magnitude values of output
%  sequence 'Windowed_Spec' can be used to accurately
%  estimate the peak amplitudes of sinusoidal components
%  in time sequence 'x'.

%   Input: 'Spec' (a sequence of complex FFT sample values)
%
%   Output: 'Windowed_Spec' (a sequence of complex flat-top
%                            windowed FFT sample values)
%
%   Based on Lyons': "Reducing FFT Scalloping Loss Errors
%   Without Multiplication", IEEE Signal Processing Magazine,
%   DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
%   Richard Lyons [December, 2011]

N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];

% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);

Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);

% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);

Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);

% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function``````

## 16QAM Modem model (Basic)

December 12, 2011 Coded in Matlab
``````clear all; close all;

n = 2^16;   %number of symbols
Fsym = 12.5; %Msps
Fs = 100;   %MHz, IF sampling frequency
Fc = 20;    %MHz, upconverter frequency

%generate 16QAM raw symbols
alphabet = [-1 -1/3 +1/3 +1];
tx = complex(randsrc(1,n,alphabet),randsrc(1,n,alphabet));

%pulse shaping, root raised cosine
h = firrcos(50,Fsym/4,.15,Fsym,'rolloff','sqrt');
tx_up1 = zeros(1,2*n);
tx_up1(1:2:end-1) = tx;
tx_shaped = filter(2*h,1,tx_up1);

%further upsampling by 4
tx_up2 = resample(tx_shaped,4,1);

%upconvert to 10MHz centre
f1 = exp(1i*2*pi*(0:8*n-1)*Fc/Fs);
tx_upconverted = tx_up2 .* f1;

%channel signal
rx_real = real(tx_upconverted);

%Rx shifts signal back to zero
f2 = exp(1i*2*pi*(0:8*n-1)*-Fc/Fs);
rx_downconverted = rx_real.*f2;

%Rx downsamples back by 4
rx_dn1 = resample(rx_downconverted,1,4);

%rx matched filter
rx_dn2 = filter(2*h,1,rx_dn1);
rx_dn2 = rx_dn2(50:end); %remove initial zeros

%one phase should be correct(odd/even)
rx_e = rx_dn2(1:2:end); %odd
rx_o = rx_dn2(2:2:end); %even

subplot(3,1,1);plot(real(tx),imag(tx),'o');
axis([-2 +2 -2 +2]);grid;
title('Tx constellations')
subplot(3,1,2);plot(real(rx_e),imag(rx_e),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, first phase')
subplot(3,1,3);plot(real(rx_o),imag(rx_o),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, second phase')``````