function Demo(L)
%Author: J. Selva.
%Date: April 2011.
%This demo constructs trials of an undamped exponential plus noise, and then computes the
%maximum likelihood estimate of its frequency, amplitude, and phase. Only one FFT of size
%2^ceil(log2(st.M)+1) is computed, plus a few operations. The complexity is that given by
%the FFT. The algoritm is explained in
%
% http://arxiv.org/abs/1104.3069
%
format long g
format compact
M = 1024; %Number of samples.
gdB = 20; %SNR at matched filter output.
ARef = 1.3; %Amplitude.
fRef = 0.372; %Frequency in [-0.5,0.5[.
PhiRef = 0.43; %Phase (rad).
if nargin > 0
for k = 1:2:length(L)
eval([L{k} '=L{k+1};']);
end
end
%Define a random vector with 1024 elements.
SNRdB = gdB - 10*log10(M); %Signal-to-noise ratio (dB).
%Complex sinusoidal amplitude, frequency, and phase.
NVar = ARef^2*10.^(-SNRdB/10); %Noise variance.
m1 = -ceil(M/2);
while 1
disp(' ')
disp(['M = ', num2str(M), ', gdB = ', num2str(gdB), ' dB, ARef = ', num2str(ARef), ...
', fRef = ', num2str(fRef), ', PhiRef = ', num2str(PhiRef), ' rad'])
disp(' ')
%Set up a trial.
a = ARef*exp(1i*2*pi*fRef*(m1:m1+M-1).'+1i*PhiRef)+(randn(M,1)+1i*randn(M,1))*sqrt(NVar/2);
%Call the nufft function. The output is a struct with some variables for later use.
st = nufft(a);
%Compute the DFT peak.
[f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st);
%The output is the following,
%
% f1, A1, and Phi1 are the ML estimates of frequency, amplitude and phase, respectively.
%
% L1 is the ML cost function value at f1 and its first two derivatives. L1(1) should be
% large since it is the square value of the peak amplitude. L1(2) is the cost
% function derivative, and it should be small since we are at a critical
% point. Finally L1(2) should be a negative number with large amplitude since f1 is
% the abcissa of a maximum (second derivative).
%
% Deg: If Deg==1 then the Newton iteration was cut too many times.
% nIt: Number of Newton iterations.
disp('Exact and estimated frequency, amplitude, and phase:')
disp([fRef,ARef,PhiRef;f1,A1,Phi1])
disp(['Number of iterations: ' num2str(nIt)])
%
%Plot the spectrum close to the ML estimate
%
f = fRef+(-10:.05:10)/st.K;
sp = nufft(st,f);
plot(f,20*log10(abs(sp)/M))
hold on
stem(f1,20*log10(A1),'r')
hold off
grid on
set(gca,'YLim',[-40,10])
xlabel('f')
ylabel('Amplitude (dB)')
text(fRef+2/st.K,5,{['Freq. ' num2str(f1)],['Ampl. ' num2str(A1)],...
['Phase. (rad) ' num2str(Phi1)]})
%
R = input('Press Enter for a new trial, q to quit: ','s');
if any(strcmp(R,{'q','Q'}))
break
end
end
%------------------------------------------------------
%------------------------------------------------------
function [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st)
%This function finds out the peak of a DFT. Its input is the struct produced by the nufft function.
%Find the spectral peak in the frequency grid.
[pr11,ind] = max(abs(st.aF).^2);
%Compute its frequency.
f0 = (ind-1)/st.K;
f0 = f0+ceil(-0.5-f0);
%L0 contains the differentials of the ML cost function of orders 0,1, 2 ...
%at f0.
L0 = LDiff(nufft(st,f0,2));
nItMax = 20;
err = 10.^-12;
nIt = 0;
%Call Iter (Newton iteration) until there is no change or the solution is degenerate.
while 1
[Deg,pr11,f1,L1] = Iter(st,f0,L0);
if Deg || abs(f1-f0) <= err || nIt >= nItMax
break
end
f0 = f1;
L0 = L1;
nIt = nIt+1;
end
pr = nufft(st,f1);
A1 = abs(pr)/st.M;
Phi1 = angle(pr);
function L = LDiff(c)
%Computes the differentials of the ML cost function from the nufft differentials.
L = [real(c(1)*conj(c(1))); 2*real(conj(c(1))*c(2)); ...
2*real(conj(c(2))*c(2)+conj(c(1))*c(3))];
function [Deg,nCut,f1,L1] = Iter(st,f0,L0)
%Performs the Newton iteration, cutting the Newton direction if necessary.
f1 = f0-L0(2)/L0(3); %Newton iteration.
f1 = f1+ceil(-0.5-f1); %Make it fall in [-0.5,0.5].
L1 = LDiff(nufft(st,f1,2)); %Compute the differentials at f1.
nCutMax = 20; %The difference f1-f0 will be reduced in the sequel if there is no ...
%improvement in the cost function
nCut = 0;
d = 1;
%In this loop the difference f1-f0 is halved a maximum of nCutMax times until there is ...
%an increase in the cost function
while L1(1) < L0(1) && nCut < nCutMax
d = d/2;
f1 = f0-d*L0(2)/L0(3);
f1 = f1+ceil(-0.5-f1);
L1 = LDiff(nufft(st,f1,2));
nCut = nCut+1;
end
Deg = nCut >= nCutMax;
%----------------------------------------------------------------
%----------------------------------------------------------------
%Non-uniform FFT function. This code is discussed in
% http://www.dsprelated.com/showcode/159.php
function varargout = nufft(varargin)
%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, see
%
%J. Selva, 'Design of barycentric interpolator for uniform and nonuniform sampling
% grids', IEEE Trans. on Signal Processing, vol. 58, n. 3, pp. 1618-1627, March 2010.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
% st = nufft(a); %a is the vector of samples.
%
% The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
% 1/st.K.
%
% The DFT is defined using
%
% A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
% where m_1=-ceil(M/2) and M the number of samples.
%
% -Then call nufft with sintax
%
% out = nufft(st,f,nd,method)
%
% where
%
% f: Vector of frequencies where the DFT is evaluated. Its elements must follow
% abs(f(k))<=0.5
%
% nd: Derivative order. nufft computes derivatives up to order nd. At the output
% out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
% column of out contains the zero-order derivatives, i.e, the values of the DFT at
% frequencies in vector f.
% method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
% +'baryLoop': The output is interpolated using the barycentric method and a loop
% implementation.
%
% +'baryVec': The output is interpolated using the barycentric method and a
% vectorized implementation.
%
% +'directLoop': The output is computed using the DFT sum directly and a loop
% implementation.
%
% +'directVec': The output is computed using the DFT sum directly and a vectorized
% implementation.
if ~isstruct(varargin{1})
st = [];
st.a = varargin{1};
st.a = st.a(:);
st.M = length(st.a);
st.m1 = -ceil(st.M/2);
st.K = 2^ceil(log2(st.M)+1);
st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
%this would reduce st.P. The number of interpolation samples is 2*st.P+1.
st.T = 1/st.K;
st.B = -2*st.m1;
st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
st.vt = MidElementToEnd((-st.P:st.P)*st.T);
st.alpha = MidElementToEnd(baryWeights(st.T,st.B,st.P));
varargout = {st};
return
end
[st,f] = deal(varargin{1:2});
nd = 0;
if nargin > 2
nd = varargin{3};
end
method = 'baryVec';
if nargin > 3
method = varargin{4};
end
Nf = length(f);
out = zeros(Nf,nd+1);
switch method
case 'baryLoop' %Interpolated method using loops
for k = 1:length(f)
x = f(k);
n = floor(x/st.T+0.5);
u = x - n * st.T;
da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');
out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);
end
case 'baryVec' %Vectorized interpolated method
f = f(:);
Nf = length(f);
n = floor(f/st.T+0.5);
u = f - n * st.T;
pr = [-st.P:-1 , 1:st.P , 0];
ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
if length(f) == 1
ms = ms.';
end
out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
case 'directLoop' % Direct method using loops
for k = 1:length(f)
out(k,1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
end
for kd = 1:nd
out(k,kd+1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,kd+1) = out(k,kd+1) + ...
((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
end
end
end
case 'directVec' %Vectorized direct method
for k = 1:length(f)
out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
for kd = 1:nd
out(k,kd+1) = ...
((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
end
end
end
varargout = {out};
function v = MidElementToEnd(v)
ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];
function v = APPulse(t,B,TSL)
v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));
function w = baryWeights(T,B,P)
vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;
N = length(vt);
LD = ones(1,N);
for k = 1:N
LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end
w = gam ./ LD;
w = w / max(abs(w));
function out = DerBarycentricInterp3(alpha,s,t,tau,nd)
vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,nd+1);
for k = 1:Nt-1
vD = vD*(tau-t(k))+alpha(k)*LF(k);
LF(k+1) = LF(k)*(tau-t(k));
end
vD = vD * (tau-t(Nt)) + alpha(Nt) * LF(Nt);
z = s;
for kd = 0:nd
z1 = z-z(end); cN = 0;
for k = 1:Nt-1
cN = cN * (tau-t(k))+z1(k)*alpha(k)*LF(k);
end
cN = cN/vD;
ztau = z(end)+(tau-t(end))*cN;
out(kd+1) = ztau;
if kd < nd
z = [ (z(1:end-1)-ztau)./(t(1:end-1)-tau) , cN ] * (kd+1);
end
end
function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)
NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
%usage.
NBlock = 1+floor(NtauL/LBtau);
Nt = length(t);
out = zeros(NtauL,nd+1);
for r = 0:NBlock-1
ind1 = 1+r*LBtau;
ind2 = min([(r+1)*LBtau,NtauL]);
Ntau = ind2-ind1+1;
z = zL(ind1:ind2,:);
tau = tauL(ind1:ind2);
vD = zeros(Ntau,1);
LF = [1,zeros(1,Nt-1)];
LF = LF(ones(Ntau,1),:);
for k = 1:Nt-1
vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
LF(:,k+1) = LF(:,k) .* (tau-t(k));
end
vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);
for kd = 0:nd
pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);
for k = 1:Nt-1
cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
end
cN = cN ./ vD;
ztau = z(:,end)+(tau-t(end)) .* cN;
out(ind1:ind2,kd+1) = ztau;
if kd < nd
pr1 = ztau;
pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
pr2 = t(1:end-1);
pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
z = [pr1 ./ pr2, cN] * (kd+1);
end
end
end
clc, clear all, close all
N=128;
x=-1:2/128:1-2/128;
%number of stages
M=q_log2(N);
%to array elements in a proper way for the butterfly
for k=0:N-1;
x1(k+1)=x(bi2de(fliplr(de2bi(k,M)))+1);
end
%to set up twiddle factors
for k=0:N/2-1;
W(k+1)=exp(-1i*2*pi*k/N);
end
%entering to graph
%k holds the stage
for k=1:M;
k2=2^k;
k21=2^(k-1);
Mk=2^(M-k);
%l holds the butterfly
for l=1:N/k2;
%l=1:N/pow(2,k)
%m holds the upper element of the butterfly
for m=1:k21;
op1 = x1((l-1)*(k2)+m);
op2 = x1((l-1)*(k2)+m+k21);
%butterfly equation
x1((l-1)*(k2)+m) = op1+W(Mk*(m-1)+1)*op2;
x1((l-1)*(k2)+m+k21) = op1-W(Mk*(m-1)+1)*op2;
end
end
end
% ****************************************************************
% baseband equivalent source of local oscillator with phase noise
% Markus Nentwig, 18.12.2011
% ****************************************************************
function pn_generator()
close all;
% ****************************************************************
% PN generator configuration
% ****************************************************************
srcPar = struct();
srcPar.n = 2 ^ 18; % generated number of output samples
srcPar.rate_Hz = 7.68e6; % sampling rate
srcPar.f_Hz = [0, 10e3, 1e6, 9e9]; % phase noise spectrum, frequencies
srcPar.g_dBc1Hz = [-80, -80, -140, -140]; % phase noise spectrum, magnitude
srcPar.spursF_Hz = [300e3, 400e3, 700e3]; % discrete spurs (set [] if not needed)
srcPar.spursG_dBc = [-50, -55, -60]; % discrete spurs, power relative to carrier
% ****************************************************************
% run PN generator
% ****************************************************************
s = PN_src(srcPar);
if false
% ****************************************************************
% export phase noise baseband-equivalent signal for use in simulator etc
% ****************************************************************
tmp = [real(s); imag(s)] .';
save('phaseNoiseSample.dat', 'tmp', '-ascii');
end
if exist('spectrumAnalyzer', 'file')
% ****************************************************************
% spectrum analyzer configuration
% ****************************************************************
SAparams = struct();
SAparams.rate_Hz = srcPar.rate_Hz; % sampling rate of the input signal
SAparams.pRef_W = 1e-3; % unity signal represents 0 dBm (1/1000 W)
SAparams.pNom_dBm = 0; % show 0 dBm as 0 dB;
SAparams.filter = 'brickwall';
SAparams.RBW_window_Hz = 1000; % convolve power spectrum with a 1k filter
SAparams.RBW_power_Hz = 1; % show power density as dBc in 1 Hz
SAparams.noisefloor_dB = -250; % don't add artificial noise
SAparams.logscale = true; % use logarithmic frequency axis
% plot nominal spectrum
figure(1); grid on;
h = semilogx(max(srcPar.f_Hz, 100) / 1e6, srcPar.g_dBc1Hz, 'k+-');
set(h, 'lineWidth', 3);
hold on;
spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'fig', 1);
ylabel('dBc in 1 Hz');
legend('nominal PSD', 'output spectrum');
title('check match with nominal PSD');
spectrumAnalyzer('signal', s, SAparams, 'fMin_Hz', 100, 'RBW_power_Hz', 'sine', 'fig', 2);
title('check carrier level (0 dBc); check spurs level(-50/-55/-60 dBc)');
ylabel('dBc for continuous-wave signal');
spectrumAnalyzer('signal', s, SAparams, 'fig', 3, 'logscale', false);
ylabel('dBc in 1 Hz');
end
end
function pn_td = PN_src(varargin)
def = {'includeCarrier', true, ...
'spursF_Hz', [], ...
'spursG_dBc', [], ...
'fMax_Hz', []};
p = vararginToStruct(def, varargin);
% length of signal in the time domain (after ifft)
len_s = p.n / p.rate_Hz
% FFT bin frequency spacing
deltaF_Hz = 1 / len_s
% construct AWGN signal in the frequency domain
% a frequency domain bin value of n gives a time domain power of 1
% for example ifft([4 0 0 0]) => 1 1 1 1
% each bin covers a frequency interval of deltaF_Hz
mag = p.n;
% scale "unity power in one bin" => "unity power per Hz":
% multiply with sqrt(deltaF_Hz):
mag = mag * sqrt(deltaF_Hz);
% Create noise according to mag in BOTH real- and imaginary value
mag = mag * sqrt(2);
% both real- and imaginary part contribute unity power => divide by sqrt(2)
pn_fd = mag / sqrt(2) * (randn(1, p.n) + 1i * randn(1, p.n));
% frequency vector corresponding to the FFT bins (0, positive, negative)
fb_Hz = FFT_freqbase(p.n, deltaF_Hz);
% interpolate phase noise spectrum on frequency vector
% note: interpolate dB on logarithmic frequency axis
H_dB = interp1(log(p.f_Hz+eps), p.g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');
% dB => magnitude
H = 10 .^ (H_dB / 20);
% H = 1e-6; % sanity check: enforce flat -120 dBc in 1 Hz
% apply filter to noise spectrum
pn_fd = pn_fd .* H;
% set spurs
for ix = 1:numel(p.spursF_Hz)
fs = p.spursF_Hz(ix);
u = abs(fb_Hz - fs);
ix2 = find(u == min(u), 1);
% random phase
rot = exp(2i*pi*rand());
% bin value of n: unity (carrier) power (0 dBc)
% scale with sqrt(2) because imaginary part will be discarded
% scale with sqrt(2) because the tone appears at both positive and negative frequencies
smag = 2 * p.n * 10 ^ (p.spursG_dBc(ix) / 20);
pn_fd(ix2) = smag * rot;
end
% limit bandwidth (tool to avoid aliasing in an application
% using the generated phase noise signal)
if ~isempty(p.fMax_Hz)
pn_fd(find(abs(fb_Hz) > p.fMax_Hz)) = 0;
end
% convert to time domain
pn_td = ifft(pn_fd);
% discard imaginary part
pn_td = real(pn_td);
% Now pn_td is a real-valued random signal with a power spectral density
% as specified in f_Hz / g_dBc1Hz.
% phase-modulate to carrier
% note: d/dx exp(x) = 1 near |x| = 1
% in other words, the phase modulation has unity gain for small phase error
pn_td = exp(i*pn_td);
if ~p.includeCarrier
% remove carrier
% returns isolated phase noise component
pn_td = pn_td - 1;
end
end
% returns a vector of frequencies corresponding to n FFT bins, when the
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
fb_Hz = 0:(n - 1);
fb_Hz = fb_Hz + floor(n / 2);
fb_Hz = mod(fb_Hz, n);
fb_Hz = fb_Hz - floor(n / 2);
fb_Hz = fb_Hz * deltaF_Hz;
end
% *************************************************************
% helper function: Parse varargin argument list
% allows calling myFunc(A, A, A, ...)
% where A is
% - key (string), value (arbitrary) => result.key = value
% - a struct => fields of A are copied to result
% - a cell array => recursive handling using above rules
% *************************************************************
function r = vararginToStruct(varargin)
% note: use of varargin implicitly packs the caller's arguments into a cell array
% that is, calling vararginToStruct('hello') results in
% varargin = {'hello'}
r = flattenCellArray(varargin, struct());
end
function r = flattenCellArray(arr, r)
ix=1;
ixMax = numel(arr);
while ix <= ixMax
e = arr{ix};
if iscell(e)
% cell array at 'key' position gets recursively flattened
% becomes struct
r = flattenCellArray(e, r);
elseif ischar(e)
% string => key.
% The following entry is a value
ix = ix + 1;
v = arr{ix};
% store key-value pair
r.(e) = v;
elseif isstruct(e)
names = fieldnames(e);
for ix2 = 1:numel(names)
k = names{ix2};
r.(k) = e.(k);
end
else
e
assert(false)
end
ix=ix+1;
end % while
end
function [time,amp] = maxima(signal, thresh, len)
% Finds the local maximum values of the given signal.
%
% Usage: [IND,AMP] = MAXIMA(X, THRESH, N);
%
% X is the data
% THRESH is the threshold that signal must be greater
% than in order to be counted. (Default = 0.0)
% N is the number of points that need to be evaluated
% (Defaults = LENGTH(X))
% IND contains the indices of X that contain the maxima data
% AMP contains the maxima amplitudes
%
% Author: sparafucile17 10/02/2003
if(exist('len') == 0)
len = length(signal);
end
if(exist('thresh') == 0)
thresh = 0.0;
end
%Initialize data
index = 1;
prev = signal(1);
local_time = 0;
local_amp = 0;
prev_slope = 1; %allow a maxima point at the second sample
%prevent length error
if(len > length(signal))
len = length(signal)
end
%Loop through data and find local maximums
for loop_i=2:len
cur = signal(loop_i);
slope = cur - prev;
if((cur < prev) & (prev_slope > 0) & (cur > thresh)) %Positive slope and amplitude dips
local_time(index) = loop_i-1;
local_amp(index) = prev;
index = index+1;
end
prev = cur;
prev_slope = slope;
end
%After loop assign data to output variables
time = local_time;
amp = local_amp;
% ********************************************
% least-mean-squares FIR design algorithm
% Markus Nentwig, 2010-2011
% release 2011/8/22
% ********************************************
function LMSFIR()
close all;
h1 = demo1('basic'); compareDemo1WithRemez(h1);
% h1 = demo1('basicLMS'); disp('demo: convergence failure is on purpose to show LMS solution');
% demo1('allpass');
% demo1('stopband');
% demo1('equalize');
% demo1('nominalResponse');
% demo1('multipassband');
% demo1('complex');
% demo1('rootRaisedCosineUpsampler');
% demo1('componentModel');
% demo1('componentModel2');
end
function h = demo1(nameOfDemo)
dpar = struct();
% parameters for a basic FIR lowpass filter design.
% kept in a struct(), so that individual examples
% can easily change them.
% sampling rate at the input of the filter
dpar.inRate_Hz = 104e6;
% number of physical FIR taps
% in a polyphase decimator, the number of internal
% coefficients will be fDecim * fStages
dpar.mStages = 36;
% end of passband. A single value will be internally
% expanded to [-9e6, 9e6].
% Asymmetric designs require
% the complexValued = true option.
% This 'default' passband can be omitted entirely, if passbands
% are declared individually later
dpar.req_passbandF_Hz = 9e6;
% defines the maximum allowed ripple in the passband.
dpar.req_passbandRipple_dB = 0.083;
% as alternative to ripple, the in-band error
% vector magnitude (EVM) can be limited
% dpar.req_passbandEVM_dB = -44;
% the passband specification may use multiple points
% dpar.req_passbandF_Hz = [-1, 5e6, 6e6, 9e6];
% dpar.req_passbandEVM_dB = [-44, -44, -34, -34];
% start of default stopband.
% as with the passband, the default stopband can be omitted,
% if individual bands are placed later.
dpar.req_stopbandF_Hz = 14e6;
dpar.req_stopbandMaxGain_dB = -30;
% dpar.req_stopbandF_Hz = [14e6, 34e6];
% dpar.req_stopbandGainMax_dB = [-30, -20];
% ********************************************
% create a filter design object "design"
% * access with LMSFIR_stage2 functions
% * evaluate with LMSFIR_stage3 function
% ********************************************
switch nameOfDemo
case 'basic'
% ********************************************
% simple filter using the parameters above
% ********************************************
design = LMSFIR_stage1_setup(dpar);
case 'basicLMS'
% ********************************************
% LMS design for comparison:
% Iterations are disabled
% ********************************************
dpar.nIter = 1;
% balance in-band / out-of-band performance as needed
dpar.inbandWeight_initValue = 5;
dpar.outOfBandWeight_initValue = 1;
design = LMSFIR_stage1_setup(dpar);
case 'allpass'
% ********************************************
% allpass design Offset the nominal delay by 1/3
% of a sample, compared to the "basic" example
% (compare the impulse responses)
% ********************************************
dpar.delayOffset = 1/3; % signal arrives now earlier
design = LMSFIR_stage1_setup(dpar);
case 'stopband'
% ********************************************
% Filter with added stopbands
% ********************************************
% the following features require more taps
dpar.mStages = 48;
% create filter design object
design = LMSFIR_stage1_setup(dpar);
% place a stopband from 14 to 16 MHz with -50 dB
design = LMSFIR_stage2_placeStopband(...
design, ...
'f_Hz', [14e6, 16e6], ...
'g_dB', [-50, -50]);
% place another stopband from 16 to 28 MHz with
% -50 dB, linearly relaxing to -40 dB
design = LMSFIR_stage2_placeStopband(...
design, ...
'f_Hz', [16e6, 28e6], ...
'g_dB', [-50, -40]);
case 'equalize'
% ********************************************
% Equalize the frequency response of another
% filter in the passband(s)
% ********************************************
% As an equalizer, this is rather inefficient with so much
% unused bandwidth. Should operate at the smallest possible BW instead.
dpar.mStages = 52;
[ffilter_Hz, H] = getExampleLaplaceDomainFilter();
% set the frequency points...
dpar.equalizeFreq_Hz = ffilter_Hz;
% ... and the filter response. The design routine will
% use linear interpolation, therefore provide a sufficiently
% dense grid.
% Equalizing an asymmetric response requires
% the complexValued=true option, leading to a complex-valued
% FIR filter.
% The equalization function needs to be normalized.
% Otherwise, pass- and stopband targets will be offset
% by the gain mismatch.
dpar.equalizeComplexGain = H;
% as alternative to the complex gain, a magnitude response
% can be given via an equalizeGain_dB argument.
% dpar.equalizeGain_dB = 20*log10(abs(H));
% an asymmetric (non-linear-phase) impulse response may
% require a group delay that is not centered in the
% FIR length.
dpar.delayOffset = 2;
design = LMSFIR_stage1_setup(dpar);
case 'componentModel'
% ********************************************
% Create a FIR filter that approximates the passband behavior of
% the analog filter accurately, and gives a similar stopband rejection
%
% The most straightforward way to model an infinite-impulse-response
% lowpass is to simply sample the impulse response. However, it needs to be
% cut to size (since the FIR filter has only finite length)
% => Chopping it off destroys the out-of-band performance (=rectangular window)
% => use a window function that trades off between passband accuracy and
% stopband rejection
% => or use the design example below.
% ********************************************
dpar.mStages = 52;
[ffilter_Hz, H] = getExampleLaplaceDomainFilter();
% set the frequency points...
dpar.nominalFreq_Hz = ffilter_Hz;
dpar.nominalComplexGain = H;
dpar.req_stopbandF_Hz = [15e6, 30e6, 55e6];
dpar.req_stopbandMaxGain_dB = [-38, -80, -115];
dpar.req_passbandF_Hz = 9e6;
% offset the impulse response, it is not centered
dpar.delayOffset = 18;
design = LMSFIR_stage1_setup(dpar);
case 'componentModel2'
% ********************************************
% an extension of "componentModel1"
% stopband rejection does not matter, but we need
% phase-accurate modeling on a region of the stopband edge
% ********************************************
dpar.mStages = 80; % this won't be cheap...
[ffilter_Hz, H] = getExampleLaplaceDomainFilter();
dpar.nominalFreq_Hz = ffilter_Hz;
dpar.nominalComplexGain = H;
dpar.req_stopbandF_Hz = [ 16e6, 100e6];
dpar.req_stopbandMaxGain_dB = [ -30, -30];
dpar.req_passbandF_Hz = 9e6;
% offset the impulse response, it is not centered
dpar.delayOffset = dpar.mStages / 2 - 8;
design = LMSFIR_stage1_setup(dpar);
% place a passband in the area on the slope that is to be modeled accurately
design = LMSFIR_stage2_placePassband(...
design, ...
'f_Hz', [12e6, 16e6], ...
'EVM_dB', [-40, -50] - 30); % nominal gain -40..-50 dB, -30 dBc EVM
case 'nominalResponse'
% ********************************************
% Design a filter to a given frequency response
% ********************************************
dpar.mStages = 50;
% the frequency response is approximated in any
% declared passband, but must be valid for any
% frequency to allow plotting.
dpar.nominalFreq_Hz = [0, 3e6, 9e6, 1e9];
dpar.nominalGain_dB = [0, 0, -6, -6];
% instead, nominalComplexGain can be used
% g = [0, 0, -3, -3];
% dpar.nominalComplexGain = 10 .^ (g/20);
design = LMSFIR_stage1_setup(dpar);
case 'multipassband'
% ********************************************
% Design a filter with three passbands
% ********************************************
dpar.mStages = 50;
dpar = rmfield(dpar, 'req_passbandF_Hz');
dpar = rmfield(dpar, 'req_passbandRipple_dB');
design = LMSFIR_stage1_setup(dpar);
design = LMSFIR_stage2_placePassband(...
design, ...
'f_Hz', [-2e6, 2e6], ...
'EVM_dB', -45);
design = LMSFIR_stage2_placePassband(...
design, ...
'f_Hz', [3e6, 7e6], ...
'EVM_dB', [-45, -40]);
design = LMSFIR_stage2_placeStopband(...
design, ...
'f_Hz', [11.8e6, 12.4e6], ...
'g_dB', -70);
case 'complex'
% ********************************************
% Design a complex-valued filter
% ********************************************
% this is also an example for what can go wrong:
% In the unconstrained section around -40 MHz, the
% frequency response is allowed to go berserk. Which
% it does.
% Solution: Place a "soft" stopband (for example at -5 dB)
% in the "don't-care" regions and add a couple of taps.
% remove passband from default parameters
dpar = rmfield(dpar, 'req_passbandF_Hz');
dpar = rmfield(dpar, 'req_passbandRipple_dB');
% remove stopband from default parameters
dpar = rmfield(dpar, 'req_stopbandF_Hz');
dpar = rmfield(dpar, 'req_stopbandMaxGain_dB');
dpar.complexValued = true;
design = LMSFIR_stage1_setup(dpar);
design = LMSFIR_stage2_placeStopband(...
design, ...
'f_Hz', [-30e6, -16e6], ...
'g_dB', -50);
design = LMSFIR_stage2_placePassband(...
design, ...
'f_Hz', [-8e6, -2e6], ...
'EVM_dB', -45);
design = LMSFIR_stage2_placeStopband(...
design, ...
'f_Hz', [3e6, 40e6], ...
'g_dB', [-30, -50]);
case 'rootRaisedCosineUpsampler'
% ********************************************
% root-raised cosine upsampling filter for WCDMA transmitter
% The input chip stream arrives at 3.84 Msps, using the
% full bandwidth.
% Before the filter, it is upsampled (zero insertion) to
% 7.68 Msps.
% The filter applies RRC-filtering with 1.22 rolloff.
% ********************************************
% calculate nominal RRC response for lookup table / linear
% interpolation
f_Hz = logspace(1, 8, 10000); f_Hz(1) = -1;
c = abs(f_Hz * 2 / 3.84e6);
c = (c-1)/(0.22); % -1..1 in the transition region
c=min(c, 1);
c=max(c, -1);
RRC_h = sqrt(1/2+cos(pi/2*(c+1))/2);
% ********************************************
% once the targets are achieved, use the remaining
% 'degrees of freedom' for least-squares optimization.
% The LMS solver will improve, where it is 'cheapest'
% The parameters are not a real-world design
% (0.5 percent EVM => -46 dB)
% ********************************************
ci = 0;
% ci = 10; % for comparison: force equiripple
dpar = struct(...
'inRate_Hz', 3.84e6, ...
'fInterp', 2, ...
'mStages', 45, ...
'req_passbandF_Hz', 3.84e6 * 1.22 / 2, ...
'req_passbandEVM_dB', -46, ...
'req_stopbandF_Hz', 2.46e6, ...
'req_stopbandMaxGain_dB', -50, ...
'nominalFreq_Hz', f_Hz, ...
'nominalGain_dB', 20*log10(RRC_h + 1e-19), ...
'convergedIterations', ci);
design = LMSFIR_stage1_setup(dpar);
[h, status] = LMSFIR_stage3_run(design);
% save('hRRC_upsampler.txt', 'h', '-ascii');
disp(status);
otherwise
assert(false);
end % switch nameOfDemo
% ********************************************
% Design the filter
% ********************************************
% h is the impulse response (FIR tap coefficients).
[h, status] = LMSFIR_stage3_run(design);
disp(status);
end
function compareDemo1WithRemez(hLMS)
% identical target settings to demo1 "basic".
% note, the demo uses targets that are exactly "on the edge"
% what the algorithm can achieve. This results in an equiripple-
% design that can be compared with remez().
% If the targets are too loosely set, pass- and stopband quality
% start to "sag" in the band middle (LMS solution => lowest overall
% error, the optimizer improves where it's "cheapest").
r_Hz = 104e6;
m = 35; % definition differs by 1
f = [0 9e6 14e6 r_Hz/2] / (r_Hz/2);
a = [1 1 0 0];
ripple_dB = 0.1;
att_dB = 30;
err1 = 1 - 10 ^ (-ripple_dB / 20);
err2 = 10 ^ (-att_dB / 20);
w = [1/err1 1/err2];
% get remez design impulse response
hRemez = remez(m, f, a, w);
figure(); hold on;
handle = plot(hLMS, 'b+'); set(handle, 'lineWidth', 3);
plot(hRemez, 'k+'); set(handle, 'lineWidth', 2);
legend('this algorithm', 'Remez');
title('comparison with Remez design (optimum equiripple)');
end
% Gets the frequency response of an "analog" (Laplace-domain) filter.
% => Chebyshev response
% => 6th order
% => 1.2 dB ripple
% => cutoff frequency at 10 MHz
% returns
% f_Hz: list of frequencies
% H: complex-valued H(f_Hz)
function [f_Hz, H] = getExampleLaplaceDomainFilter()
[IIR_b, IIR_a] = cheby1(6, 1.2, 1, 's');
% evaluate it on a wide enough frequency range
f_Hz = logspace(1, 10, 1000); f_Hz(1) = -1;
% Laplace domain operator for normalized frequency
fCutoff_Hz = 10e6;
s = 1i * f_Hz / fCutoff_Hz;
% polynomial in s
H_num = polyval(IIR_b, s);
H_denom = polyval(IIR_a, s);
H = H_num ./ H_denom;
end
% === LMSFIR_xyz "API" functions ===
% ********************************************
% LMSFIR_stagex_... functions to interact with design 'object'
% to be executed in 'stage'-order
% ********************************************
function d = LMSFIR_stage1_setup(varargin)
p = varargin2struct(varargin);
d = struct();
% number of frequency points. Increase to improve accuracy.
% Frequencies are quantized to +/- rate / (2 * nSamples)
d.nSamples = 1024;
% default polyphase interpolation: none
d.fInterp = 1;
% default polyphase decimation: none
d.fDecim = 1;
% max. number of iterations
d.nIter = 100;
% for pure LMS solution, set nIter to 1 and change the weights below as needed
d.inbandWeight_initValue = 1;
d.outOfBandWeight_initValue = 1;
% abort when the iteration weights grow too large.
% This happens when targets are impossible.
% The result may still be meaningful, though.
d.abortWeight = 1e12;
% keep iterating, if the targets are reached.
% Once the "equi"-ripple iteration has brought all peaks to an acceptable level,
% the LMS solver will use the remaining "degrees of freedom" for a LMS optimization.
% The solver improves "where it is easy / cheap". This results in sloped
% stopbands and "drooping" EVM in passbands.
% Often, LMS is the best choice => set converged iterations to 0.
d.convergedIterations = 10;
% for a complex-valued filter, use "true".
% With a real-valued design, user input is only evaluated for positive frequencies!
d.complexValued = false;
% by default, the basis waveforms given to the optimizer are
% within a delay range of +/- half the FIR length.
% For nonlinear phase types (equalization / nominal frequency
% response), this may be suboptimal.
% Meaningful values shouldn't exceed +/- half the number of
% coefficients in the impulse response.
d.delayOffset = 0;
% copy parameters
fn = fieldnames(p);
for ix = 1:size(fn, 1)
key = fn{ix};
d.(key) = p.(key);
end
% frequency base over FFT range
fb = 0:(d.nSamples - 1);
fb = fb + floor(d.nSamples / 2);
fb = mod(fb, d.nSamples);
fb = fb - floor(d.nSamples / 2);
fb = fb / d.nSamples; % now [0..0.5[, [-0.5..0[
fb = fb * d.inRate_Hz * d.fInterp;
d.fb = fb;
% in real-valued mode, negative frequencies are treated as
% positive, when 'user input' is evaluated
if d.complexValued
d.fbUser = fb;
else
d.fbUser = abs(fb);
end
% ********************************************
% target settings. Those will be modified by
% LMSFIR_stage2_xyz()
% ********************************************
% initial value of NaN indicates: all entries are unset
d.errorSpecBinVal_inband_dB = zeros(size(d.fb)) + NaN;
d.errorSpecBinVal_outOfBand_dB = zeros(size(d.fb)) + NaN;
% ********************************************
% process req_passband requirement
% needs to be done at stage 1, because it is
% used for 'gating' with the tightenExisting /
% relaxExisting options
% ********************************************
if isfield(d, 'req_passbandF_Hz')
par = struct('onOverlap', 'error');
if isfield(d, 'req_passbandRipple_dB')
par.ripple_dB = d.req_passbandRipple_dB;
end
if isfield(d, 'req_passbandEVM_dB')
par.EVM_dB = d.req_passbandEVM_dB;
end
par.f_Hz = d.req_passbandF_Hz;
d = LMSFIR_stage2_placePassband(d, par);
end % if req_passbandF_Hz
% ********************************************
% process req_stopband requirement
% needs to be done at stage 1, because it is
% used for 'gating' with the tightenExisting /
% relaxExisting options
% ********************************************
if isfield(d, 'req_stopbandF_Hz')
f_Hz = d.req_stopbandF_Hz;
g_dB = d.req_stopbandMaxGain_dB;
% extend to infinity
if isscalar(f_Hz)
f_Hz = [f_Hz 9e19];
g_dB = [g_dB g_dB(end)];
end
d = placeBand...
(d, ...
'f_Hz', f_Hz, 'g_dB', g_dB, ...
'type', 'stopband', ...
'onOverlap', 'tighten');
end
% ********************************************
% plot management
% ********************************************
d.nextPlotIx = 700;
end
function d = LMSFIR_stage2_placeStopband(d, varargin)
p = varargin2struct(varargin);
% shorthand notation g_dB = -30; f_Hz = 9e6;
% extend fixed passband to positive infinity
if isscalar(p.f_Hz)
assert(p.f_Hz > 0);
p.f_Hz = [p.f_Hz 9e99];
end
if isscalar(p.g_dB)
p.g_dB = ones(size(p.f_Hz)) * p.g_dB;
end
% default action is to use the stricter requirement
if ~isfield(p, 'onOverlap')
p.onOverlap = 'tighten';
end
d = placeBand(d, 'type', 'stopband', p);
end
function d = LMSFIR_stage2_placePassband(d, varargin)
p = varargin2struct(varargin);
% default action is to use the stricter requirement
if ~isfield(p, 'onOverlap')
p.onOverlap = 'tighten';
end
% translate ripple spec to error
if isfield(p, 'ripple_dB')
assert(p.ripple_dB > 0);
eSamplescale = 10 ^ (p.ripple_dB / 20) - 1;
EVM_dB = 20*log10(eSamplescale);
end
if isfield(p, 'EVM_dB')
EVM_dB = p.EVM_dB;
end
% convert scalar to two-element vector
if isscalar(EVM_dB)
EVM_dB = [EVM_dB EVM_dB];
end
% *** handle f_Hz ***
f_Hz = p.f_Hz;
% extend to 0 Hz
if isscalar(f_Hz)
f_Hz = [0 f_Hz];
end
% *** create the passband ***
d = placeBand(d, ...
'type', 'passband', ...
'f_Hz', f_Hz, ...
'g_dB', EVM_dB, ...
'onOverlap', p.onOverlap);
end
% ********************************************
% the filter design algorithm
% h: impulse response
% status: converged or not
% note that even if convergence was not reached,
% the resulting impulse response is "the best we
% can do" and often meaningful.
% ********************************************
function [h, status] = LMSFIR_stage3_run(d)
1;
% mTaps is number of physical FIR stages
% m is number of polyphase coefficients
d.m = d.mStages * d.fInterp;
% masks flagging pass-/stopband frequencies
mask_inband = NaN_to_0_else_1(d.errorSpecBinVal_inband_dB);
mask_outOfBand = NaN_to_0_else_1(d.errorSpecBinVal_outOfBand_dB);
% sanity check... (orthogonality of wanted and unwanted component)
assert(sum(mask_inband) > 0, 'passband is empty');
assert(sum(mask_inband .* mask_outOfBand) == 0, ...
'passband and stopband overlap');
% ********************************************
% start with flat passband signals at input and output of filter
% those will become the input to the LMS solver.
% ********************************************
sigSolverAtInput_fd = mask_inband;
sigSolverAtOutput_fd = sigSolverAtInput_fd;
% ********************************************
% for even-sized FFT length, there is one bin at the
% Nyquist limit that gives a [-1, 1, -1, 1] time domain
% waveform. It has no counterpart with opposite frequency
% sign and is therefore problematic (time domain component
% cannot be delayed).
% Don't assign any input power here.
% ********************************************
if mod(d.nSamples, 2) == 0
ixNyquistBin = floor(d.nSamples/2) + 1;
sigSolverAtInput_fd(ixNyquistBin) = 0;
sigSolverAtOutput_fd(ixNyquistBin) = 0;
end
if isfield(d, 'equalizeFreq_Hz')
% ********************************************
% Filter equalizes a given passband frequency response
% ********************************************
if isfield(d, 'equalizeGain_dB')
cgain = 10 .^ (equalizeGain_dB / 20);
else
cgain = d.equalizeComplexGain;
end
d.Heq = evaluateFilter(d.fb, d.equalizeFreq_Hz, cgain, d.complexValued);
assert(isempty(find(isnan(d.Heq), 1)), ...
['equalizer frequency response interpolation failed. ' ...
'Please provide full range data for plotting, even if it does not ', ...
'affect the design']);
% ********************************************
% apply frequency response to input signal.
% The LMS solver will invert this response
% ********************************************
sigSolverAtInput_fd = sigSolverAtInput_fd .* d.Heq;
end
if isfield(d, 'nominalFreq_Hz')
% ********************************************
% (equalized) filter matches a given passband frequency response
% ********************************************
if isfield(d, 'nominalGain_dB')
cgain = 10 .^ (d.nominalGain_dB / 20);
else
cgain = d.nominalComplexGain;
end
d.Hnom = evaluateFilter(d.fb, d.nominalFreq_Hz, cgain, d.complexValued);
assert(isempty(find(isnan(d.Hnom), 1)), ...
['nominal frequency response interpolation failed. ' ...
'Please provide full range data for plotting, even if it does not ', ...
'affect the design']);
% ********************************************
% apply frequency response to output signal.
% The LMS solver will adapt this response
% ********************************************
sigSolverAtOutput_fd = sigSolverAtOutput_fd .* d.Hnom;
end
% ********************************************
% compensate constant group delay from equalizer and nominal
% frequency response. This isn't optimal, but it is usually
% a good starting point (use delayOffset parameter)
% ********************************************
[coeff, ref_shiftedAndScaled, deltaN] = fitSignal_FFT(...
ifft(sigSolverAtInput_fd), ifft(sigSolverAtOutput_fd));
% the above function also scales for best fit. This is not desired here, instead
% let the LMS solver match the gain. Simply scale it back:
ref_shifted = ref_shiftedAndScaled / coeff;
sigSolverAtOutput_fd = fft(ref_shifted);
if false
% ********************************************
% plot time domain waveforms (debug)
% ********************************************
figure(76); hold on;
plot(fftshift(abs(ifft(sigSolverAtOutput_fd))), 'k');
plot(fftshift(abs(ifft(sigSolverAtInput_fd))), 'b');
title('time domain signals');
legend('reference (shifted)', 'input signal');
end
% ********************************************
% main loop of the design algorithm
% => initialize weights
% => loop
% => design optimum LMS filter that transforms weighted input
% into weighted output
% => adapt weights
% => iterate
% ********************************************
% at this stage, the input to the algorithm is as follows:
% => errorSpec for in-band and out-of-band frequencies
% (masks are redundant, can be derived from above)
% => LMS_in_fd and
% => LMS_out_fd: Signals that are given to the LMS solver.
% Its task is: "design a FIR filter that transforms LMS_in_fd into LMS_out_fd".
% initialize weights
outOfBandWeight = mask_outOfBand * d.outOfBandWeight_initValue;
inbandWeight = mask_inband * d.inbandWeight_initValue;
status = '? invalid ?';
hConv = [];
remConvIter = d.convergedIterations;
for iter=1:d.nIter
% inband weight is applied equally to both sides to shape the error
% out-of-band weight is applied to the unwanted signal
LMS_in_fd = sigSolverAtInput_fd .* inbandWeight...
+ mask_outOfBand .* outOfBandWeight;
LMS_out_fd = sigSolverAtOutput_fd .* inbandWeight;
% ********************************************
% cyclic time domain waveforms from complex spectrum
% ********************************************
LMS_in_td = ifft(LMS_in_fd);
LMS_out_td = ifft(LMS_out_fd);
% ********************************************
% construct FIR basis (output per coeffient)
% time domain waveforms, shifted according to each FIR tap
% ********************************************
basis = zeros(d.m, d.nSamples);
% introduce group delay target
ix1 = -d.m/2+0.5 + d.delayOffset;
ix2 = ix1 + d.m - 1;
rowIx = 1;
for ix = ix1:ix2 % index 1 appears at ix1
basis(rowIx, :) = FFT_delay(LMS_in_td, ix);
rowIx = rowIx + 1;
end
if d.complexValued
rightHandSide_td = LMS_out_td;
else
% use real part only
basis = [real(basis)];
rightHandSide_td = [real(LMS_out_td)];
pRp = real(rightHandSide_td) * real(rightHandSide_td)' + eps;
pIp = imag(rightHandSide_td) * imag(rightHandSide_td)';
assert(pIp / pRp < 1e-16, ...
['got an imaginary part where there should be none. ', ...
'uncomment the following lines, if needed']);
% if designing a real-valued equalizer for a complex-valued frequency response,
% use the following to solve LMS over the average:
% basis = [real(basis) imag(basis)];
% rightHandSide_td = [real(LMS_out_td), imag(LMS_out_td)];
end
% ********************************************
% LMS solver
% find a set of coefficients that scale the
% waveforms in "basis", so that their sum matches
% "rightHandSide_td" LMS-optimally
% ********************************************
pbasis = pinv(basis .');
h = transpose(pbasis * rightHandSide_td .');
% pad impulse response to n
irIter = [h, zeros(1, d.nSamples-d.m)];
% undo the nominal group delay
irIter = FFT_delay(irIter, ix1);
HIter = fft(irIter);
% ********************************************
% filter test signal
% ********************************************
eq_fd = sigSolverAtInput_fd .* HIter;
% ********************************************
% subtract actual output from targeted output
% results in error spectrum
% ********************************************
err_fd = sigSolverAtOutput_fd - eq_fd;
err_fd = err_fd .* mask_inband; % only in-band matters
EVM_dB = 20*log10(abs(err_fd)+1e-15);
% ********************************************
% out-of-band leakage
% ********************************************
leakage_dB = 20*log10(abs(HIter .* mask_outOfBand + 1e-15));
% ********************************************
% compare achieved and targeted performance
% ********************************************
deltaLeakage_dB = leakage_dB - d.errorSpecBinVal_outOfBand_dB;
deltaEVM_dB = EVM_dB - d.errorSpecBinVal_inband_dB;
% ********************************************
% find bins where performance should be improved
% or relaxed
% ********************************************
ixImprLeakage = find(deltaLeakage_dB > 0);
ixImprEVM = find(deltaEVM_dB > 0);
ixRelLeakage = find(deltaLeakage_dB < -3);
ixRelEVM = find(deltaEVM_dB < -3);
status = 'iteration limit reached';
if isempty(ixImprLeakage) && isempty(ixImprEVM)
% both targets met. Convergence!
if remConvIter > 0
remConvIter = remConvIter - 1;
status = 'converged once, now trying to improve';
hConv = h;
else
status = 'converged';
break;
end
end
% ********************************************
% improve / relax in-band and out-of-band
% ********************************************
if ~isempty(ixImprLeakage)
% tighten out-of-band
outOfBandWeight(ixImprLeakage) = outOfBandWeight(ixImprLeakage)...
.* 10 .^ ((deltaLeakage_dB(ixImprLeakage) + 0.1) / 20);
end
if ~isempty(ixRelLeakage)
% relax out-of-band
outOfBandWeight(ixRelLeakage) = outOfBandWeight(ixRelLeakage)...
.* 10 .^ (deltaLeakage_dB(ixRelLeakage) / 3 / 20);
end
if ~isempty(ixImprEVM)
% tighten in-band
inbandWeight(ixImprEVM) = inbandWeight(ixImprEVM)...
.* 10 .^ ((deltaEVM_dB(ixImprEVM) + 0.01) / 20);
end
if ~isempty(ixRelEVM)
% relax in-band
inbandWeight(ixRelEVM) = inbandWeight(ixRelEVM)...
.* 10 .^ (deltaEVM_dB(ixRelEVM) / 2 / 20);
end
if max([inbandWeight, outOfBandWeight] > d.abortWeight)
status = 'weight vector is diverging';
break;
end
end % for iter
% ********************************************
% recover from convergence failure after convergence
% during improvement phase
% ********************************************
if ~strcmp(status, 'converged')
if ~isempty(hConv)
h = hConv;
status = 'converged';
end
end
if true
% ********************************************
% plot impulse response
% ********************************************
if d.complexValued
figure(); hold on;
stem(real(h), 'k');
stem(imag(h), 'b');
legend('real(h)', 'imag(h)');
else
figure(); hold on;
stem(h);
legend('h');
end
title('impulse response');
end
% ********************************************
% plot frequency response
% ********************************************
inbandBins = find(mask_inband);
outOfBandBins = find(mask_outOfBand);
d=doPlotStart(d, ['Frequency response (Status:', status, ')']);
d=doPlotH(d, HIter, 'b', '|H_{design}(f)|', 2);
handle = plot(d.fb(outOfBandBins), d.errorSpecBinVal_outOfBand_dB(outOfBandBins), 'b+');
set(handle, 'markersize', 2);
d=addLegend(d, 'req. stopband');
d = doPlot_dB(d, EVM_dB, 'r', 'error');
handle = plot(d.fb(inbandBins), d.errorSpecBinVal_inband_dB(inbandBins), 'r+');
set(handle, 'markersize', 2);
d=addLegend(d, 'req. passband error');
d=doPlotEnd(d);
ylim([-100, 10]);
if false
% ********************************************
% plot constraint signal and weights
% ********************************************
figure(31); grid on; hold on;
handle = plot(fftshift(d.fb), fftshift(20*log10(mask_outOfBand)));
set(handle, 'lineWidth', 3);
x = d.fb; y = 20*log10(inbandWeight / max(inbandWeight));
handle = plot(x(inbandBins), y(inbandBins), 'k+'); set(handle, 'lineWidth', 3);
x = d.fb;
y = 20*log10(outOfBandWeight / max(outOfBandWeight));
handle = plot(x(outOfBandBins), y(outOfBandBins), 'b+'); set(handle, 'lineWidth', 3);
xlabel('f/Hz'); ylabel('dB');
ylim([-80, 40]);
legend('constraint signal', 'in-band weight', 'out-of-band weight');
title('weighting factor (normalized to 0 dB)');
end
hasEq = isfield(d, 'Heq');
hasNom = isfield(d, 'Hnom');
if hasEq || hasNom
% ********************************************
% plot equalization / nominal target
% ********************************************
d=doPlotStart(d, 'equalization / nominal target');
d=doPlotH(d, HIter, 'b', '|H_{design}(f)|', 2);
if hasEq
d=doPlotH(d, d.Heq, 'k', '|H_{eq}(f)| to equalize (invert)');
eqR = HIter .* d.Heq;
d=doPlotH(d, eqR, 'c', '|H_{design}(f)H_{eq}(f)|', 2);
handle = plot(d.fb(inbandBins), ...
20*log10(abs(eqR(inbandBins)) + 1e-15), 'c*');
set(handle, 'markersize', 3);
d=addLegend(d, '|H_{eq}(in-band)');
end
if hasNom
d = doPlotH(d, d.Hnom, 'g', '|H_{nom}|', 2);
handle = plot(d.fb(inbandBins), ...
20*log10(abs(HIter(inbandBins)) + 1e-15), 'b*');
set(handle, 'markersize', 3);
d=addLegend(d, '|H_{design}(f)H_{eq}(f) in-band');
end
d=doPlotEnd(d);
% set y-range
ymax = 20*log10(max(abs(HIter)));
ylim([-50, ymax+3]);
end
end
% === LMSFIR helper functions ===
% evaluates frequency response f_dB; g_Hz at fb
% the return value will contain NaN for out-of-range entries
% in fb
function binVal = buildBinVal(varargin)
p = varargin2struct(varargin);
f_Hz = p.f_Hz;
g_dB = p.g_dB;
% shorthand notation f = [f1, f2]; g = -30;
if isscalar(g_dB)
g_dB = ones(size(f_Hz)) * g_dB;
end
% tolerate sloppy two-argument definition
if size(f_Hz, 2) == 2 && f_Hz(1) > f_Hz(2)
f_Hz = fliplr(f_Hz);
g_dB = fliplr(g_dB);
end
binVal = interp1(f_Hz, g_dB, p.fbUser, 'linear');
end
function d = placeBand(d, varargin)
p = varargin2struct(varargin);
% create requirements vector
binVal = buildBinVal('f_Hz', p.f_Hz, ...
'g_dB', p.g_dB, ...
'fbUser', d.fbUser);
% look up requirements vector from design object
switch p.type
case 'passband'
fn = 'errorSpecBinVal_inband_dB';
case 'stopband'
fn = 'errorSpecBinVal_outOfBand_dB';
otherwise
assert(false);
end
designObject_binVal = d.(fn);
% check overlap
if strcmp(p.onOverlap, 'error')
m1 = NaN_to_0_else_1(designObject_binVal);
m2 = NaN_to_0_else_1(binVal);
assert(isempty(find(m1 .* m2, 1)), ...
['newly declared band overlaps existing band, '...
'which was explicitly forbidden by onOverlap=error']);
p.onOverlap = 'tighten'; % there won't be overlap,
% merging is dummy operation
end
% merging rules
switch p.onOverlap
case 'tighten'
logicOp = 'or';
valueOp = 'min';
case 'relax'
logicOp = 'or';
valueOp = 'max';
case 'tightenExisting'
logicOp = 'and';
valueOp = 'min';
case 'relaxExisting'
logicOp = 'and';
valueOp = 'max';
otherwise
assert(false);
end
% merge requirements tables
binValMerged = mergeBinVal(...
'binVal1', designObject_binVal, ...
'binVal2', binVal, ...
'logicalOperator', logicOp, ...
'valueOperator', valueOp);
% assign new requirements table
d.(fn) = binValMerged;
end
function r = NaN_to_0_else_1(vec)
r = zeros(size(vec));
% logical indexing, instead of r(find(~isnan(vec))) = 1;
r(~isnan(vec)) = 1;
end
function binVal = mergeBinVal(varargin)
p = varargin2struct(varargin);
% region where first argument is defined
mask1 = NaN_to_0_else_1(p.binVal1);
% region where second argument is defined
mask2 = NaN_to_0_else_1(p.binVal2);
% region where result will be defined
switch(p.logicalOperator)
case 'or'
mask = mask1 + mask2;
case 'and'
mask = mask1 .* mask2;
otherwise
assert(false);
end
ix = find(mask);
% merge into result
binVal = zeros(size(p.binVal1)) + NaN;
switch(p.valueOperator)
case 'min'
% note: The function min/max ignore NaNs (see "min" man
% page in Matlab)
% if one entry is NaN, the other entry will be returned
binVal(ix) = min(p.binVal1(ix), p.binVal2(ix));
case 'max'
binVal(ix) = max(p.binVal1(ix), p.binVal2(ix));
otherwise
assert(false);
end
end
% evaluates [f / gain] filter specification on the frequency grid
function H = evaluateFilter(f_eval, modelF, modelH, complexValued)
oneSided = false;
if ~complexValued
oneSided = true;
else
if min(modelF) > min(f_eval)
disp(['Warning: Filter model does not contain (enough) negative frequencies. ', ...
'assuming symmetric H(f) / real-valued h(t)']);
oneSided = true;
end
end
if oneSided
f_evalOrig = f_eval;
f_eval = abs(f_eval);
end
H = interp1(modelF, modelH, f_eval, 'linear');
if oneSided
% enforce symmetry (=> real-valued impulse response)
logicalIndex = (f_evalOrig < 0);
H(logicalIndex) = conj(H(logicalIndex));
end
end
function [d, handle] = doPlotH(d, H, spec, legEntry, linewidth)
handle = plot(fftshift(d.fb), fftshift(20*log10(abs(H)+1e-15)), spec);
d = addLegend(d, legEntry);
if exist('linewidth', 'var')
set(handle, 'lineWidth', linewidth);
end
end
function [d, handle] = doPlot_dB(d, H, spec, legEntry, linewidth)
handle = plot(fftshift(d.fb), fftshift(H), spec);
d.legList{size(d.legList, 2) + 1} = legEntry;
if exist('linewidth', 'var')
set(handle, 'lineWidth', linewidth);
end
end
function d = doPlotStart(d, plotTitle)
figure(d.nextPlotIx);
title(plotTitle);
grid on; hold on;
d.nextPlotIx = d.nextPlotIx + 1;
d.legList = {};
end
function d = doPlotEnd(d)
legend(d.legList);
xlabel('f/Hz');
ylabel('dB');
end
function d = addLegend(d, legEntry)
d.legList{size(d.legList, 2) + 1} = legEntry;
end
% === general-purpose library functions ===
% handling of function arguments
% someFun('one', 1, 'two', 2, 'three', 3) => struct('one', 1, 'two', 2, 'three', 3)
% a struct() may appear in place of a key ('one') and gets merged into the output.
function r = varargin2struct(arg)
assert(iscell(arg));
switch(size(arg, 2))
case 0 % varargin was empty
r=struct();
case 1 % single argument, wrapped by varargin into a cell list
r=arg{1}; % unwrap
assert(isstruct(r));
otherwise
r=struct();
% iterate through cell elements
ix=1;
ixMax=size(arg, 2);
while ix <= ixMax
e=arg{ix};
if ischar(e)
% string => key/value. The next field is a value
ix = ix + 1;
v = arg{ix};
r.(e) = v;
elseif isstruct(e)
names = fieldnames(e);
assert(size(names, 2)==1); % column
for ix2 = 1:size(names, 1)
k = names{ix2};
v = e.(k);
r.(k) = v;
end
else
disp('invalid token in vararg handling. Expecting key or struct. Got:');
disp(e);
assert(false)
end
ix=ix+1;
end % while
end % switch
end
function sig = FFT_delay(sig, nDelay)
sig = fft(sig); % to frequency domain
nSigSamples = size(sig, 2);
binFreq=(mod(((0:nSigSamples-1)+floor(nSigSamples/2)), nSigSamples)-floor(nSigSamples/2));
phase = -2*pi*nDelay / nSigSamples .* binFreq;
rot = exp(1i*phase);
if mod(nSigSamples, 2)==0
% even length - bin at Nyquist limit
rot(nSigSamples/2+1)=cos(phase(nSigSamples/2+1));
end
sig = sig .* rot;
sig = ifft(sig); % to time domain
end
% *******************************************************
% delay-matching between two signals (complex/real-valued)
%
% * matches the continuous-time equivalent waveforms
% of the signal vectors (reconstruction at Nyquist limit =>
% ideal lowpass filter)
% * Signals are considered cyclic. Use arbitrary-length
% zero-padding to turn a one-shot signal into a cyclic one.
%
% * output:
% => coeff: complex scaling factor that scales 'ref' into 'signal'
% => delay 'deltaN' in units of samples (subsample resolution)
% apply both to minimize the least-square residual
% => 'shiftedRef': a shifted and scaled version of 'ref' that
% matches 'signal'
% => (signal - shiftedRef) gives the residual (vector error)
%
% Example application
% - with a full-duplex soundcard, transmit an arbitrary cyclic test signal 'ref'
% - record 'signal' at the same time
% - extract one arbitrary cycle
% - run fitSignal
% - deltaN gives the delay between both with subsample precision
% - 'shiftedRef' is the reference signal fractionally resampled
% and scaled to optimally match 'signal'
% - to resample 'signal' instead, exchange the input arguments
% *******************************************************
function [coeff, shiftedRef, deltaN] = fitSignal_FFT(signal, ref)
n=length(signal);
% xyz_FD: Frequency Domain
% xyz_TD: Time Domain
% all references to 'time' and 'frequency' are for illustration only
forceReal = isreal(signal) && isreal(ref);
% *******************************************************
% Calculate the frequency that corresponds to each FFT bin
% [-0.5..0.5[
% *******************************************************
binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
% *******************************************************
% Delay calculation starts:
% Convert to frequency domain...
% *******************************************************
sig_FD = fft(signal);
ref_FD = fft(ref, n);
% *******************************************************
% ... calculate crosscorrelation between
% signal and reference...
% *******************************************************
u=sig_FD .* conj(ref_FD);
if mod(n, 2) == 0
% for an even sized FFT the center bin represents a signal
% [-1 1 -1 1 ...] (subject to interpretation). It cannot be delayed.
% The frequency component is therefore excluded from the calculation.
u(length(u)/2+1)=0;
end
Xcor=abs(ifft(u));
% figure(); plot(abs(Xcor));
% *******************************************************
% Each bin in Xcor corresponds to a given delay in samples.
% The bin with the highest absolute value corresponds to
% the delay where maximum correlation occurs.
% *******************************************************
integerDelay = find(Xcor==max(Xcor));
% (1): in case there are several bitwise identical peaks, use the first one
% Minus one: Delay 0 appears in bin 1
integerDelay=integerDelay(1)-1;
% Fourier transform of a pulse shifted by one sample
rotN = exp(2i*pi*integerDelay .* binFreq);
uDelayPhase = -2*pi*binFreq;
% *******************************************************
% Since the signal was multiplied with the conjugate of the
% reference, the phase is rotated back to 0 degrees in case
% of no delay. Delay appears as linear increase in phase, but
% it has discontinuities.
% Use the known phase (with +/- 1/2 sample accuracy) to
% rotate back the phase. This removes the discontinuities.
% *******************************************************
% figure(); plot(angle(u)); title('phase before rotation');
u=u .* rotN;
% figure(); plot(angle(u)); title('phase after rotation');
% *******************************************************
% Obtain the delay using linear least mean squares fit
% The phase is weighted according to the amplitude.
% This suppresses the error caused by frequencies with
% little power, that may have radically different phase.
% *******************************************************
weight = abs(u);
constRotPhase = 1 .* weight;
uDelayPhase = uDelayPhase .* weight;
ang = angle(u) .* weight;
r = [constRotPhase; uDelayPhase] .' \ ang.'; %linear mean square
%rotPhase=r(1); % constant phase rotation, not used.
% the same will be obtained via the phase of 'coeff' further down
fractionalDelay=r(2);
% *******************************************************
% Finally, the total delay is the sum of integer part and
% fractional part.
% *******************************************************
deltaN = integerDelay + fractionalDelay;
% *******************************************************
% provide shifted and scaled 'ref' signal
% *******************************************************
% this is effectively time-convolution with a unit pulse shifted by deltaN
rotN = exp(-2i*pi*deltaN .* binFreq);
ref_FD = ref_FD .* rotN;
shiftedRef = ifft(ref_FD);
% *******************************************************
% Again, crosscorrelation with the now time-aligned signal
% *******************************************************
coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
shiftedRef=shiftedRef * coeff;
if forceReal
shiftedRef = real(shiftedRef);
end
end
% *********************************************
% Weighting factors for Remez equiripple FIR design
% M. Nentwig
% *********************************************
close all; clear all;
f = [0:9]/10; % normalized frequency 0..1
% a = nominal |H(f)| on a linear scale (sample amplitude)
% 1 1 : passband from frequency 0 to 0.1
% 0 0 : stopband from frequency 0.2 to 0.3
% 0.5 0.5 : passband with -6 dB (0.5 amplitude) from frequency 0.4 to 0.5
% 0 0 : stopband from frequency 0.6 to 0.7
% 1 1 : passband from frequency 0.8 to 0.9
% other frequency ranges are "don't care" areas.
a = [1 1 0 0 0.5 0.5 0 0 1 1];
% *********************************************
% design specification for each band
% *********************************************
ripple1_dB = 0.3;
att2_dB = 60;
ripple3_dB = 0.2;
att4_dB = 70;
ripple5_dB = 0.1;
% *********************************************
% error for each band, on a linear scale
% *********************************************
% note: passband 3 has a nominal gain of 0.5 in 'a'.
% the allowed error of |H(f)| scales accordingly.
err1 = 1 - 10 ^ (-ripple1_dB / 20);
err2 = 10 ^ (-att2_dB / 20);
err3 = (1 - 10 ^ (-ripple3_dB / 20)) * 0.5;
err4 = 10 ^ (-att4_dB / 20);
err5 = 1 - 10 ^ (-ripple5_dB / 20);
% the weight of each band is chosen as the inverse of the targeted error
% (stricter design target => higher weight).
% we could normalize each entry in w to (err1+err2+err3+err4+err5)
% but it would appear as common factor in all entries and therefore make no difference.
w = [1/err1 1/err2 1/err3 1/err4 1/err5];
% filter order
% this design problem is 'tweaked' so that the resulting filter is (quite) exactly on target
% for the given n, which can be changed only in steps of 1.
% Note that increasing order by 1 can make the filter worse due to even / odd
% number of points (impulse response symmetry)
n = 52;
% *********************************************
% Run Remez / Parks McClellan filter design
% *********************************************
h = remez(n, f, a, w);
% force row vector (OctaveForge and Matlab's "remez" differ)
h = reshape(h, 1, prod(size(h)));
% *********************************************
% Plot the results
% *********************************************
figure(1); hold on;
% zero-pad the impulse response to set frequency resolution
% of the FFT
h = [h, zeros(1, 1000)];
% frequency base
npts = size(h, 2);
fbase = (0:npts-1)/npts;
plot(fbase, 20*log10(abs((fft(h)+1e-15))), 'b');
xlim([0, 0.5]);
% *********************************************
% sketch the requirements
% *********************************************
e = [1 1];
plot([f(1), f(2)]/2, e * ripple1_dB, 'k');
plot([f(1), f(2)]/2, e * -ripple1_dB, 'k');
plot([f(3), f(4)]/2, e * -att2_dB, 'k');
plot([f(5), f(6)]/2, e * ripple3_dB - 6.02, 'k');
plot([f(5), f(6)]/2, e * -ripple3_dB - 6.02, 'k');
plot([f(7), f(8)]/2, e * -att4_dB, 'k');
plot([f(9), f(10)]/2, e * ripple5_dB, 'k');
plot([f(9), f(10)]/2, e * -ripple5_dB, 'k');
xlabel('normalized frequency 0..1');
ylabel('dB');
% **************************************************************
% bank-switched Farrow resampler
% M. Nentwig, 2011
% Note: Uses cyclic signals (wraps around)
% **************************************************************
close all; clear all;
% inData contains input to the resampling process (instead of function arguments)
inData = struct();
% **************************************************************
% example coefficients.
% Each column [c0; c1; c2; ...] describes a polynomial for one tap coefficent in fractional time ft [0, 1]:
% tapCoeff = c0 + c1 * ft + c2 * ft ^ 2 + ...
% Each column corresponds to one tap.
% the matrix size may be changed arbitrarily.
%
% The example filter is based on a 6th order Chebyshev Laplace-domain prototype.
% **************************************************************
if false
% for comparison, this is a conventional design (no bank switching)
inData.cMatrix =[ -8.57738278e-3 7.82989032e-1 7.19303539e+000 6.90955718e+000 -2.62377450e+000 -6.85327127e-1 1.44681608e+000 -8.79147907e-1 7.82633997e-2 1.91318985e-1 -1.88573400e-1 6.91790782e-2 3.07723786e-3 -6.74800912e-3
2.32448021e-1 2.52624309e+000 7.67543936e+000 -8.83951796e+000 -5.49838636e+000 6.07298348e+000 -2.16053205e+000 -7.59142947e-1 1.41269409e+000 -8.17735712e-1 1.98119464e-1 9.15904145e-2 -9.18092030e-2 2.74136108e-2
-1.14183319e+000 6.86126458e+000 -6.86015957e+000 -6.35135894e+000 1.10745051e+001 -3.34847578e+000 -2.22405694e+000 3.14374725e+000 -1.68249886e+000 2.54083065e-1 3.22275037e-1 -3.04794927e-1 1.29393976e-1 -3.32026332e-2
1.67363115e+000 -2.93090391e+000 -1.13549165e+000 5.65274939e+000 -3.60291782e+000 -6.20715544e-1 2.06619782e+000 -1.42159644e+000 3.75075865e-1 1.88433333e-1 -2.64135123e-1 1.47117661e-1 -4.71871047e-2 1.24921920e-2] / 231.46 * 20;
inData.nBanks = 1;
else
% same example filter as above, but now the matrix contains three alternative coefficient banks for each tap.
% The order was reduced from cubic to quadratic.
% column 1: first bank, tap 1
% column 2: second bank, tap 1
% column 3: third bank, tap 1
% column 4: first bank, tap 2
% and so on
inData.cMatrix =[ 2.87810386e-4 4.70096244e-3 7.93412570e-2 4.39824536e-1 1.31192924e+000 2.67892232e+000 4.16465421e+000 5.16499621e+000 5.15592605e+000 3.99000369e+000 2.00785470e+000 -7.42377060e-2 -1.52569354e+000 -1.94402804e+000 -1.40915797e+000 -3.86484652e-1 5.44712939e-1 9.77559688e-1 8.32191447e-1 3.22691788e-1 -2.13133045e-1 -5.08501962e-1 -4.82928807e-1 -2.36313854e-1 4.76034568e-2 2.16891966e-1 2.20894063e-1 1.08361553e-1 -2.63421832e-2 -1.06276015e-1 -1.07491548e-1 -5.53793711e-2 4.86314061e-3 3.94357182e-2 4.06217506e-2 2.17199064e-2 1.60318761e-3 -8.40370106e-3 -8.10525279e-3 -3.62112499e-3 -4.13413072e-4 2.33101911e-4
-3.26760325e-3 -6.46028234e-3 1.46793247e-1 5.90235537e-1 1.18931309e+000 1.57853546e+000 1.40402774e+000 5.76506323e-1 -6.33522788e-1 -1.74564700e+000 -2.24153717e+000 -1.91309453e+000 -9.55568978e-1 1.58239169e-1 9.36193787e-1 1.10969783e+000 7.33284446e-1 1.06542194e-1 -4.15412084e-1 -6.06616434e-1 -4.54898908e-1 -1.20841199e-1 1.82941623e-1 3.12543429e-1 2.49935829e-1 8.05376898e-2 -7.83213666e-2 -1.47769751e-1 -1.18735248e-1 -3.70656555e-2 3.72608374e-2 6.71425397e-2 5.17812605e-2 1.55564930e-2 -1.40896327e-2 -2.35058137e-2 -1.59635057e-2 -3.44701792e-3 4.14108065e-3 4.56234829e-3 1.59503132e-3 -3.17301882e-4
5.64310141e-3 7.74786707e-2 2.11791763e-1 2.84703201e-1 1.85158633e-1 -8.41118142e-2 -3.98497442e-1 -5.86821615e-1 -5.40397941e-1 -2.47558080e-1 1.50864737e-1 4.59312895e-1 5.41539400e-1 3.84673917e-1 9.39576331e-2 -1.74932542e-1 -3.01635463e-1 -2.56239225e-1 -9.87146864e-2 6.82216764e-2 1.59795852e-1 1.48668245e-1 6.62563431e-2 -2.71234898e-2 -8.07045577e-2 -7.76841351e-2 -3.55333136e-2 1.23206602e-2 3.88535040e-2 3.64199073e-2 1.54608563e-2 -6.59814558e-3 -1.72735099e-2 -1.46307777e-2 -5.04363288e-3 3.31049461e-3 6.01267607e-3 3.83904192e-3 3.92549958e-4 -1.36315264e-3 -9.76017430e-4 7.46699178e-5] / 133.64 * 20;
inData.nBanks = 3;
end
% **************************************************************
% Create example signal
% **************************************************************
nIn = 50; % used for test signal generation only
if false
% complex signal
inData.signal = cos(2*pi*(0:(nIn-1)) / nIn) + cos(2*2*pi*(0:(nIn-1)) / nIn) + 1.5;
else
% impulse response
inData.signal = zeros(1, nIn); inData.signal(1) = 1; %
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 1) = 1; % enable to show constant c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(2, 1) = 1; % enable to show linear c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(3, 1) = 1; % enable to show quadratic c in first tap, first bank
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 2) = 1; % enable to show constant c in first tap, second bank
end
% **************************************************************
% Resample to the following number of output samples
% must be integer, otherwise arbitrary
% **************************************************************
inData.nSamplesOut = floor(nIn * 3 * 6.28);
% **************************************************************
% Set up Farrow resampling
% **************************************************************
nSamplesIn = size(inData.signal, 2);
nSamplesOut = inData.nSamplesOut;
order = size(inData.cMatrix, 1) - 1; % polynomial order
% number of input samples that contribute to one output sample (FIR size)
nTaps = size(inData.cMatrix, 2);
assert(mod(nTaps, inData.nBanks) == 0);
nTaps = nTaps / inData.nBanks; % only one out of nBanks coefficients contributes at any time
% pointer to the position in the input stream for each output sample (row vector, real numbers), starting at 0
inputIndex = (0:nSamplesOut-1) / nSamplesOut * nSamplesIn;
% split into integer part (0..nSamplesIn - 1) ...
inputIndexIntegerPart = floor(inputIndex);
% ... and fractional part [0, 1[
inputIndexFractionalPart = inputIndex - inputIndexIntegerPart;
% bank switching
% the fractional part is again split into an integer and fractional part
inputIndexFractionalPart = inputIndexFractionalPart * inData.nBanks;
inputIndexFractionalPart_int = floor(inputIndexFractionalPart); % coefficient bank index
inputIndexFractionalPart_frac = inputIndexFractionalPart - inputIndexFractionalPart_int; % fractional time 0..1 within each bank
% **************************************************************
% Calculate output stream
% First constant term (conventional FIR), then linear, quadratic, cubic, ...
% **************************************************************
outStream = zeros(1, inData.nSamplesOut);
for ixOrder = 0 : order
% note: fractional time is now defined for each sub-segment [0, 1[
x = inputIndexFractionalPart_frac .^ ixOrder;
% **************************************************************
% Add the contribution of each tap one-by-one
% **************************************************************
for ixTap = 0 : nTaps - 1
% coefficient bank switching: There are inData.nBanks alternative coefficients for each tap
c = inData.cMatrix(ixOrder+1, ixTap * inData.nBanks + inputIndexFractionalPart_int + 1);
% index of input sample that contributes to output via the current tap
% higher tap index => longer delay => older input sample => smaller data index
dataIx = inputIndexIntegerPart - ixTap;
% wrap around
dataIx = mod(dataIx, nSamplesIn);
% array indexing starts at 1
dataIx = dataIx + 1;
delayed = inData.signal(dataIx);
% for each individual output sample (index in row vector),
% - evaluate f = c(order, tapindex) * fracPart .^ order
% - scale the delayed input sample with f
% - accumulate the contribution of all taps
% this implementation performs the operation for all output samples in parallel (row vector)
outStream = outStream + c .* delayed .* x;
end % for ixTap
end % for ixOrder
% **************************************************************
% plot
% **************************************************************
xIn = linspace(0, 1, nSamplesIn + 1); xIn = xIn(1:end-1);
xOut = linspace(0, 1, nSamplesOut + 1); xOut = xOut(1:end-1);
figure(); grid on; hold on;
stem(xIn, inData.signal, 'k+-');
plot(xOut, outStream, 'b+-');
legend('input', 'output');
title('bank-switched Farrow resampling. Signals are cyclic.');
function [y] = signal_unmute(x, index, duration, Fs)
% Return the input vector with a unmute that occurs at a specific
% index and with an exponential ramp-up to reduce "pop" sounds.
% The output will start off at -100dB gain and end at 0dB gain.
%
% Usage: y = SIGNAL_UNMUTE(x, index, duration, Fs);
%
% X is your one-dimensional input array
% INDEX is where in the input signal you want the unmute to begin
% DURATION is how long (in seconds) to exponentially ramp up the
% input signal. 100ms is recommended.
% FS is the sample rate
%
% Example:
% You want to unmute your signal at 0.5sec with a duration of 100ms
% and with a 48k Fs. (signal is unmuted completely at 0.6sec)
% y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003
% Input must have some length
if(length(x) == 1)
error('ERROR: input signal must have more than one element');
end
% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
error('ERROR: Input must be one-dimensional');
end
% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
error(['There are not enough samples in X to complete the unmute. ' ...
'Either change the mute duration or move the index back.' ]);
end
% Flip vector (temporarily)
if((size(x, 2) ~= 1))
x = x';
flip_back = true;
else
flip_back = false;
end
% Calculate exponential coefficient
dB_atten = -100; %What do we consider "muted" to be in decibels
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));
% Build the Gain array
gain = [zeros(index, 1); ones((length(x)-index), 1);];
b = [ coeff 0];
a = [ 1 (-1*(1-coeff))];
gain = filter(b,a,gain);
% Apply Mute (gain) to the input signal
y = gain .* x;
% Flip the vector (if required)
if(flip_back == true);
y = y';
end
% *************************************************************
% Peak-to-average analyzis of complex baseband-equivalent signal
% Markus Nentwig, 10/12/2011
% Determines the magnitude relative to the RMS average
% which is exceeded with a given (small) probability
% A typical probability value is 99.9 %, which is very loosely related
% to an uncoded bit error rate target of 0.1 %.
% *************************************************************
function sn_headroom()
% number of samples for example signals
n = 1e6;
% *************************************************************
% example: 99.9% PAR for white Gaussian noise
% *************************************************************
signal = 12345 * (randn(1, n) + 1i * randn(1, n));
[PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, false);
printf('white Gaussian noise: %1.1f dB PAR at radio frequency\n', PAR_dB);
[PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, true);
printf('white Gaussian noise: %1.1f dB PAR at baseband\n', PAR_dB);
% *************************************************************
% example: PAR of continuous-wave signal
% *************************************************************
% a quarter cycle gives the best coverage of amplitude values
% the statistics of the remaining three quarters are identical (symmetry)
signal = 12345 * exp(1i*2*pi*(0:n-1) / n / 4);
% Note: peaks can be below the average, depending on the given probability
% a possible peak-to-average ratio < 1 (< 0 dB) is a consequence of the
% peak definition and not unphysical.
% For a continuous-wave signal, the average value equals the peak value,
% and PAR results slightly below 0 dB are to be expected.
[PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, false);
printf('continuous-wave: %1.1f dB PAR at radio frequency\n', PAR_dB);
[PAR, PAR_dB] = peakToAverageRatio(signal, 0.999, true);
printf('continuous-wave: %1.1f dB PAR at baseband\n', PAR_dB);
% *************************************************************
% self test:
% For a real-valued Gaussian random variable, the probability
% to not exceed n sigma is
% n = 1: 68.27 percent
% n = 2: 95.5 percent
% n = 3: 99.73 percent
% see http://en.wikipedia.org/wiki/Normal_distribution
% section "confidence intervals"
% *************************************************************
signal = 12345 * randn(1, n);
[PAR, PAR_dB] = peakToAverageRatio(signal, 68.2689/100, false);
printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(1), PAR_dB);
[PAR, PAR_dB] = peakToAverageRatio(signal, 95.44997/100, false);
printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(2), PAR_dB);
[PAR, PAR_dB] = peakToAverageRatio(signal, 99.7300/100, false);
printf('expecting %1.3f dB PAR, got %1.3f dB\n', 20*log10(3), PAR_dB);
end
function [PAR, PAR_dB] = peakToAverageRatio(signal, peakProbability, independentIQ)
1;
% force row vector
assert(min(size(signal)) == 1, 'matrix not allowed');
signal = signal(:) .';
assert(peakProbability > 0);
assert(peakProbability <= 1);
% *************************************************************
% determine RMS average and normalize signal to unity
% *************************************************************
nSamples = numel(signal);
sAverage = sqrt(signal * signal' / nSamples);
signal = signal / sAverage;
% *************************************************************
% "Peaks" in a complex-valued signal can be defined in two
% different ways:
% *************************************************************
if ~independentIQ
% *************************************************************
% Radio-frequency definition:
% ---------------------------
% The baseband-equivalent signal represents the modulation on
% a radio frequency carrier
% The instantaneous magnitude of the modulated radio frequency
% wave causes clipping, for example in a radio frequency
% power amplifier.
% The -combined- magnitude of in-phase and quadrature signal
% (that is, real- and imaginary part) is relevant.
% This is the default definition, when working with radio
% frequency (or intermediate frequency) signals, as long as a
% single, real-valued waveform is being processed.
% *************************************************************
sMag = abs(signal);
t = 'mag(s) := |s|; cdf(mag(s))';
else
% *************************************************************
% Baseband definition
% -------------------
% The baseband-equivalent signal is explicitly separated into
% an in-phase and a quadrature branch that are processed
% independently.
% The -individual- magnitudes of in-phase and quadrature branch
% cause clipping.
% For example, the definition applies when a complex-valued
% baseband signal is processed in a digital filter with real-
% valued coefficients, which is implemented as two independent,
% real-valued filters on real part and imaginary part.
% *************************************************************
% for each sample, use the maximum of in-phase / quadrature.
% If both clip in the same sample, it's counted only once.
sMag = max(abs(real(signal)), abs(imag(signal)));
t = 'mag(s) := max(|real(s)|, |imag(s)|); cdf(mag(s))';
end
% *************************************************************
% determine number of samples with magnitude below the "peak"
% level, based on the given peak probability
% for example: probability = 0.5 => nBelowPeakLevel = nSamples/2
% typically, this will be a floating point number
% *************************************************************
nBelowPeakLevel = peakProbability * nSamples;
% *************************************************************
% interpolate the magnitude between the two closest samples
% *************************************************************
sMagSorted = sort(sMag);
x = [0 1:nSamples nSamples+1];
y = [0 sMagSorted sMagSorted(end)];
magAtPeakLevel = interp1(x, y, nBelowPeakLevel, 'linear');
% *************************************************************
% Peak-to-average ratio
% signal is normalized, average is 1
% the ratio relates to the sample magnitude
% "power" is square of magnitude => use 20*log10 for dB conversion.
% *************************************************************
PAR = magAtPeakLevel;
PAR_dB = 20*log10(PAR + eps);
% *************************************************************
% illustrate the CDF and the result
% *************************************************************
if true
figure();
plot(y, x / max(x));
title(t);
xlabel('normalized magnitude m');
ylabel('prob(mag(s) < m)');
line([1, 1] * magAtPeakLevel, [0, 1]);
line([0, max(y)], [1, 1] * peakProbability);
end
end
% ------------------------------------------------------
% Title: formfilter.m
%
% Author: David Valencia
%
% Institution: UPIITA-IPN 2010
%
% for DSPrelated.com
%
% Published in: http://www.dsprelated.com/showcode/12.php
%
% Description: This program receives 2 basic ortogonal
% filters (one is high-pass and one is low-pass)
% the filtering level and the branch number.
% As a result it returns the equivalent filter
% of the specified branch.
%
% Dependencies: upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revision: v1.0a
% - Commented and translated in English
%
% For more details on this code and its implementation
% see the following blog posts:
%
% http://www.dsprelated.com/showarticle/115.php
% http://www.dsprelated.com/showarticle/116.php
%
% ------------------------------------------------------
function [hx] = formfilter(n_stages,branch,h0,h1)
p = branch;
% Seed vector
hx = 0;
hx(1) = 1;
switch n_stages
case 1
% If it is a single stage filter
% the branches are just the base filters
if mod(branch,2) ~= 0
hx = h0;
else
hx = h1;
end
case 2
% For a 2 stage filter bank, the
% equivalent filters are simply
% the convolution between the corresponding
% base filters, and one of them is upsampled
% The use of upsample2 is needed to avoid having
% a large quantitiy of zeros at the end of the vector
% which certainly difficults its usage to build a
% convolution matrix.
switch branch
case 1
hx = conv(h0,upsample2(h0,2));
case 2
hx = conv(h0,upsample2(h1,2));
case 3
hx = conv(h1,upsample2(h0,2));
case 4
hx = conv(h1,upsample2(h1,2));
otherwise
beep;
fprintf('\nFor a 2 stage filter bank there can not be a fifth branch');
end
otherwise
% For a n>2 stages filter bank, a more ellaborated
% process must be made. A series of upsamplings and convolutions
% are needed to get an equivalent vector.
for i=0:n_stages-2
q = floor(p /(2^(n_stages-1-i)));
if (q == 1)
hx = conv(hx,upsample2(h1,2^i));
else
hx = conv(hx,upsample2(h0,2^i));
end
p = mod(p,2^(n_stages-1-i));
end
% Depending on the parity of the branch number, the filter
% goes through a last convolution
t = mod(branch,2);
if(t == 1)
hx = conv(hx,upsample2(h0,2^(n_stages-1)));
else
hx = conv(hx,upsample2(h1,2^(n_stages-1)));
end
end
% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% by means of a linear convolution matrix
%
% Author: David Valencia
% UPIITA IPN
%
% Posted at: http://www.dsprelated.com/showcode/11.php
%
% Description:
% Computes the Discrete Wavelet Transform
% of n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% convolution matrix method (linear convolution)
%
% Dependencies:
%
% formfilter.m
% http://www.dsprelated.com/showcode/12.php
% upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revisions:
% v1.0a Commented and translated in English
%
% ----------------------------------------------------
close all; clear all; clc;
disp('== GENERALIZED MATRICIAL DWT ==')
%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db3';
% Obtain base filters
if(typeofbasis == 'b')
[Rf,Df] = biorwavf(typbior);
[h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
[h0,h1,g0,g1] = wfilters(typor);
end;
%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.
N = 16; %Number of input samples
x = (1:N)';
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
dec_fact = n_branches; %Decimating factor
% L = Basic filters lenght (h0 รณ h1)
% N = Input vector length
% n_stages = number of stages, it generates (2^n_stages) branches
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
fprintf('\nLinear convolution matrix processing\n');
%% STEP 3: Build analysis matrices
% -------------------------------------------------------------
% High-pass analysis matrix (single scale) [linear convolution]
% As the DWT has a simple high-pass branch, the matrix is easy
% to build with the base filter.
% Append zeros to each side of the vector to use it as a sliding
% part to build the matrix
movil_1 = [zeros(1,N-1),fliplr(h1),zeros(1,N-1)];
lm = length(movil_1);
n_rows = ceil((L+N-1)/2); %Compute the number of rows
W1 = zeros(n_rows,N); % Allocate memory for the high pass matrix
dec_factor = 2; %Decimate factor for the high-pass matrix
% Build matrix (with sequential shifting)
for i = 0:(n_rows-1)
W1(i+1,:) = movil_1(lm-(i*dec_factor)-N+1 : lm-(i*dec_factor)); %Rama pasaaltas
end
disp('High-pass analysis matrix: ');
W1
disp('High-pass branch results');
y1 = W1 * x
% -------------------------------------------------------------
% Low-pass several stage analysis matrix [linear convolution]
% Get an equivalent filter to get its length
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);
% Compute the number of rows per branch
rowsperfilter = ceil((N+Lhx)/dec_fact);
% Allocate memory for the low-pass matrix
W0 = zeros((n_branches/2)*rowsperfilter,N);
tic; %Start cronometer
% Build low pass filter
for i = 0:(n_branches/2)-1
% Get the equivalent filter depending on the number of stages
if n_stages < 3
hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
else
hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
end
% Append zeros to the vector
movil = [zeros(1,N-1),fliplr(hx),zeros(1,N-1)];
lm = length(movil);
% Shift and add the vector to the matrix
for j = 0:rowsperfilter-1
W0(i*rowsperfilter + j + 1,:) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
end
end
disp('Low-pass many stage analysis matrix: ');
W0
disp('Low-pass filter bank results: ');
y0 = W0*x
%% STEP 4: Recover signal
disp('Recovered signal: ');
xn1 = (W1')*y1;
xn0 = (W0')*y0;
xn = xn1+xn0
%% STEP 5: Compute error
err = (x - xn).^2;
disp('Error :');
err = (sum(err)/N);
disp(err);
fprintf('\nProcessing time = %d seg \n',toc);
fprintf('\nProgram end >>> \n');
% David Valencia
% Posted @ DSPRelated.com
% in: http://www.dsprelated.com/showcode/10.php
function [y] = upsample2(x,N)
lx = length(x);
y = 0;
% Put first sample in the first position of the output
y(1) = x(1);
for i=2:lx
% Create a "chunk" (or piece) of vector, composed
% with a series of zeros and the next sample from
% the input
chunk = [zeros(1,N-1),x(i)];
% Append the chunk to the current output vector
y = [y,chunk];
end
function xdft=fdft(x)
N=length(x);
xdft(1,N)=0; % defining a 1xN sequence with all the elements equals to zero.
%*************THIS IS THE DFT IMPLEMENTATION**************
%For each value of k, the Summatory function is computed over the interval 0<n<N-1.
for k=0:N-1,
for n=0:N-1,
xdft(k+1)=xdft(k+1) + x(n+1)*exp(-j*2*pi*k*n/N);
end
end
%NOTE: You may skip the use of for-loops by defining k and n as two N-points sequences. This is:
% k=0:N-1; n=0:N-1;
%The DFT is now computed in matrix form (as a matrix multiplication of the N-points Fourier Matrix and the discrete-time sequence x).
% *****************************************************************
% Using the Moore-Penrose Pseudoinverse to remove known disturbances from
% measured data
%
% all vectors are column vectors.
% *****************************************************************
close all; clear all;
% *****************************************************************
% Example: A measured signal has been polluted with power-line hum,
% at a fundamental frequency, 3rd and 5th harmonic, with unknown
% amplitude and phase.
% *****************************************************************
m=1000; % number of measured points
fs=22000; % sampling rate
fHum=50; % known 50 Hz power line hum frequency
index=(1:m)';
% create power-line hum terms
hum=0.87654*sin(2*pi*index*fHum/fs+0.12345);
hum3rd=0.65432*sin(2*pi*index*3*fHum/fs+0.45678);
hum5th=0.3321*sin(2*pi*index*5*fHum/fs+0.78901);
% create actual data
idealData=0.1*randn(m, 1);
% create measured data - deteriorated by hum
measuredData=idealData+hum+hum3rd+hum5th;
% *****************************************************************
% Processing to remove the hum
% *****************************************************************
% We know frequency and number of harmonics
% We DO NOT the amplitude and phase of each term.
% note: a*sin(ometaT)+b*cos(omegaT) is equivalent to c*sin(omegaT+d)
% where either [a, b] or [c, d] are the unknown parameters.
basis=[sin(2*pi*index*1*fHum/fs), cos(2*pi*index*1*fHum/fs), ...
sin(2*pi*index*3*fHum/fs), cos(2*pi*index*3*fHum/fs), ...
sin(2*pi*index*5*fHum/fs), cos(2*pi*index*5*fHum/fs)];
% *****************************************************************
% Moore-Penrose Pseudoinverse: Least-squares fit between the basis
% waveforms and the measured signal.
% *****************************************************************
leastSquaresCoeff=pinv(basis)*measuredData;
reconstructedHum=basis*leastSquaresCoeff;
correctedData=measuredData-reconstructedHum;
% *****************************************************************
% Plot measured signal and reconstructed hum
% *****************************************************************
figure(); hold on;
plot(measuredData);
plot(reconstructedHum, 'r', 'LineWidth', 2);
legend({'measured data', 'reconstructed hum signal'});
% *****************************************************************
% Plot the remainder after subtracting the reconstructed hum.
% Compare with ideal result
% *****************************************************************
figure(); hold on;
plot(idealData, 'k', 'LineWidth', 2);
plot(correctedData);
plot(idealData-correctedData, 'r', 'LineWidth', 3);
legend({'ideal result', 'measured and corrected', 'error'})
% *****************************************************************
% Audio interpolation using Cubic Hermite Spline (Catmull-Rom)
%
% The audio file samples describe a continuous-waveform.
% The task is to interpolate between samples, for example in sample rate
% conversion or variable-speed playback.
% The cubic spline approach is one particular trade-off between
% quality and efficiency.
%
% For an efficient C-implementation, refer to fluid_rvoice_dsp.c
% in the fluidsynth project on sourceforge.
% Based on work by Olli Niemitalo:
% http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf page 11
% ******************************************************************
% build coefficient matrix
% ******************************************************************
% ideally, the interpolation coefficients would be recalculated for any possible
% value of ptrFrac below.
% example input file: http://www.elisanet.fi/mnentwig/webroot/common/in.wav
% Note: This example is written for clarity, not speed.
nPoly=128;
x=((0:nPoly-1)/nPoly)'; % column vector
c1=(x .* (-0.5 + x .* (1 - 0.5 .* x)));
c2=(1.0 + x .* x .* (1.5 .* x - 2.5));
c3=(x .* (0.5 + x .* (2.0 - 1.5 .* x)));
c4=(0.5 .* x .* x .* (x - 1.0));
% ******************************************************************
% load input data
% ******************************************************************
[y, fs, bits]=wavread('in.wav');
[nSamples, nChan]=size(y);
inputdata=transpose(y(:, 1)); % first channel if stereo
% prepare output data storage
outIndex=1; outputData=[]; nAlloc=0;
% ******************************************************************
% "ptr" indicates a position in the input data.
% Loop through the whole range until end-of-data.
% ******************************************************************
ptr=2;
while ptr < nSamples-2
% "rate" determines the number of input samples per output sample
% rate=1.1; % faster constant-speed playback
% variable playback speed (arbitrary example)
rate=1+0.3*sin(outIndex/40000);;
% ******************************************************************
% Interpolate samples
% ******************************************************************
% ptr is the index into the input data, needed for the next output sample.
% Usually, it lies betwen points, for example ptr=1234.56
% integer part (1234)
ptrBase=floor(ptr);
% fractional part [0..1[, for example 0.56
ptrFrac=ptr-ptrBase;
% look up the filter for the fractional part
fIndex=floor(ptrFrac*nPoly)+1;
% input samples
s1=inputdata(ptrBase-1);
s2=inputdata(ptrBase);
s3=inputdata(ptrBase+1);
s4=inputdata(ptrBase+2);
% performance optimization. Reallocate some extra space.
if outIndex>=nAlloc
nAlloc=nAlloc+2000;
outputData(nAlloc)=0;
end
% ******************************************************************
% Calculate output sample.
% ******************************************************************
outputData(outIndex)=s1*c1(fIndex)+s2*c2(fIndex)+s3*c3(fIndex)+s4*c4(fIndex);
% that's it. Next output sample
outIndex=outIndex+1;
% advance the index in the source data
ptr=ptr+rate;
end
% performance optimization part II, chop extra memory
outputData=outputData(1:outIndex-1)';
% ******************************************************************
% Write output to file
% ******************************************************************
% wavwrite('out.wav', outputData, fs, bits); % OCTAVE version
wavwrite(outputData, fs, bits, 'out.wav'); % MATLAB version
% *************************************************************
% Raised cosine window for OFDM intersymbol transitions
% *************************************************************
% Use the window for the guard period between adjacent symbols
% to smoothen the waveform transition, limit sidelobes and prevent
% unwanted emissions into adjacent frequency ranges
%
% reference:
% Data transmission by Frequency Division Multiplexing
% using the Discrete Fourier Transform
% S.B. Weinstein and Paul M. Ebert
% IEEE transactions on Communications technology, No. 5,
% October 1971
% *************************************************************
% Length of the guard period in samples
n=64;
index=linspace(0, 1, n+2);
index=index(2:end-1);
windowNextSymbol=(1-cos(pi*index))/2;
windowPreviousSymbol=fliplr(windowNextSymbol);
% plot
figure(); title('raised cosine window for OFDM symbol transitions');
xlabel('sample index into the guard period between symbols');
ylabel('window weight');
hold on;
plot(windowPreviousSymbol, 'r+');
plot(windowNextSymbol, 'b+');
legend({'cyclic extension of symbol i-1', 'cyclic extension of symbol i'});
% ******************************************************
% Symbol error rate of Quadrature Amplitude Modulation
% ******************************************************
% uncoded, additive white Gaussian noise channel
% Implements [1] Figure 5.2-16 on a wider SNR range
% [1] John G. Proakis
% "Digital Communications"
% 4th edition, McGraw-Hill", 2001
% http://www.amazon.com/Digital-Communications-John-Proakis/dp/0072321113
% Note that different definitions are in common use:
% a) {symbol, bit} errors
% b) signal-to-noise ratio per {symbol, bit}
close all; clear all;
SNR_perSymbol_dB=[-3:0.1:42];
% definition: Below [1] 5.2-78
SNR_perSymbol=10.^(SNR_perSymbol_dB/10);
plotPerBitX=[]; plotY=[];
% number of constellation points (M-ary QAM)
for M=[4, 16, 64, 256, 1024]
% the energy of each symbol is used to transmit log2(M) bits
SNR_perBit=SNR_perSymbol/log2(M);
% [1] 5.2-78 argument of Q(...)
Qarg=sqrt(3/(M-1)*SNR_perSymbol);
% [1] 2.1-98
Q=1/2*erfc(Qarg/sqrt(2));
% [1] eq. 5.2-77
% probability of error for PAM per quadrature signal
PsqrtM=2*(1-1/sqrt(M))*Q;
% [1] eq. 5.2-79
PM=1-(1-PsqrtM).^2;
plotPerBitX=[plotPerBitX; 10*log10(SNR_perBit)];
plotY=[plotY; PM];
end
figure(1);
h = semilogy(SNR_perSymbol_dB', plotY'); grid on;
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per symbol (dB)'); ylabel('symbol error rate');
xlim([0, 40]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});
figure(2);
h = semilogy(plotPerBitX', plotY'); grid on; ylabel('symbol error rate');
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per bit (dB)');
xlim([-3, 30]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});