## Phase drift of NCO

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

n = 50000;                            %number of test samples
Fs = 100;                             %sampling rate in Msps
Fo = 1.17;                            %target NCO frequency in MHz
M = [2^16 10000];                     %NCO phase resolution, two cases
for i = 1:2
inc = round(M(i)*Fo/Fs);                         %phase increment
lut = round(2047*cos(2*pi*(0:M(i)-1)/M(i)));     %lut, one cycle cosine data
addr = 0:inc:inc*(n-1);
addr = round(mod(addr,M(i)));
addr(addr >= M(i)) = addr(addr >= M(i)) - M(i);  %check address overflow
y(i,1:n) = lut(addr+1);                          %add 1 for matlab LUT
end

fprintf('increment value = %2.6f\r',M*Fo/Fs);

subplot(2,1,1);hold;
plot([y(1,1:1000) zeros(1,100) y(1,end-999:end)],'.-');
plot([y(2,1:1000) zeros(1,100) y(2,end-999:end)],'r.--');
legend('drift','no drift');
title('initial and last 1000 samples');

[P1, F1] = pwelch(y(1,:), hann(n), 0, n, Fs);
[P2, F2] = pwelch(y(2,:), hann(n), 0, n, Fs);
subplot(2,1,2); hold
plot(F1,10*log10(P1));
plot(F2,10*log10(P2),'r--');
legend('drift','no drift')
axis ([0 50 -20 100]);
grid
xlabel('MHz')``````

## Spectrum analyzer

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

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

noise = randn(1, n);
pulse = zeros(1, n); pulse(1) = 1;

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

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

% ************************************
% example (logarithmic frequency scale,
% no RBW filter)
% ************************************
fPar = struct('rate_Hz', rate_Hz, ...
'chebyOrder', 6, ...
'chebyRipple_dB', 1, ...
'fCorner_Hz', 1e5);

saPar = struct('rate_Hz', rate_Hz, ...
'pRef_W', 1e-3, ...
'fMin_Hz', 1e3, ...
'RBW_power_Hz', 1e5, ...
'RBW_window_Hz', 1e3, ...
'filter', 'none', ...
'plotStyle', 'b-', ...
'logscale', true, ...
'fig', 2, ...
'noisefloor_dB', -300);

tmp1 = normalize(IIR_filterExample('signal', noise, fPar));
tmp2 = normalize(IIR_filterExample('signal', pulse, fPar));

spectrumAnalyzer('signal', tmp1, saPar);
spectrumAnalyzer('signal', tmp2, saPar, 'plotStyle', 'k-');
end

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

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

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

n = numel(p.signal);
fb_Hz = FFT_frequencyBasis(n, p.rate_Hz);

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

% clip to -1..1 range
c=min(c, 1);
c=max(c, -1);

% raised-cosine filter
H = 1/2+cos(pi/2*(c+1))/2;

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

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

% ************************************
% continuous-time filter example
% ************************************
function sig = IIR_filterExample(varargin)
p = vararginToStruct(varargin);

% design continuous-time IIR filter
[IIR_b, IIR_a] = cheby1(p.chebyOrder, p.chebyRipple_dB, 1, 's');

% evaluate on the FFT frequency grid
fb_Hz = FFT_frequencyBasis(numel(p.signal), p.rate_Hz);

% Laplace domain operator, normalized to filter cutoff frequency
% (the above cheby1 call designs for unity cutoff)
s = 1i * fb_Hz / p.fCorner_Hz;

% polynomial in s
H_num = polyval(IIR_b, s);
H_denom = polyval(IIR_a, s);
H = H_num ./ H_denom;

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

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

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

handle=nan; % avoid warning

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

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

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

scale_to_dBm=sqrt(p.pRef_W/0.001);

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

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

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

winLen=0;
if ~isempty(winLeft)

% zero-pad the power spectrum in the middle with a length
% equivalent to the window size.
% this guarantees that there is enough bandwidth for the filter!
% one window length is enough, the spillover from both sides overlaps
% Scale accordingly.
winLen=size(winLeft, 2)+size(winRight, 2);

% note: fr is the power shown in the plot, NOT a frequency
% domain representation of a signal.
% There is no need to renormalize because of the length change
center=floor(nSamples/2)+1;
rem=nSamples-center;
fr=[fr(1:center-1), zeros(1, winLen-1), fr(center:end)];
% construct window with same length as fr
win=zeros(size(fr));
win(1:1+size(winLeft, 2)-1)=winLeft;
win(end-size(winRight, 2)+1:end)=winRight;

assert(isreal(win));
assert(isreal(fr));
assert(size(win, 2)==size(fr, 2));

% convolve using FFT
win=fft(win);
fr=fft(fr);

fr=fr .* win;
fr=ifft(fr);
fr=real(fr); % chop off roundoff error imaginary part
fr=max(fr, 0); % chop off small negative numbers

% remove padding
fr=[fr(1:center-1), fr(end-rem+1:end)];
end

% ************************************
% build frequency basis and rotate 0 Hz to center
% ************************************
fb = FFT_frequencyBasis(numel(fr), p.rate_Hz);
fr = fftshift(fr);
fb = fftshift(fb);

if false
% use in special cases (very long signals)

% ************************************
% data reduction:
% If a filter is used, details smaller than windowBW are
% suppressed already by the filter, and using more samples
% gives no additional information.
% ************************************
if numel(fr) > p.NMax
switch(p.filter)
case 'gaussian'
% 0.2 offset from the peak causes about 0.12 dB error
incr=floor(windowBW/binspacing_Hz/5);
case 'brickwall'
% there is no error at all for a peak
incr=floor(windowBW/binspacing_Hz/3);
case 'none'
% there is no filter, we cannot discard data at this stage
incr=-1;
end

if incr > 1
fr=fr(1:incr:end);
fb=fb(1:incr:end);
end
end
end

% ************************************
% data reduction:
% discard beyond fMin / fMax, if given
% ************************************
indexMin = 1;
indexMax = numel(fb);
flag=0;
if isfield(p, 'fMin_Hz')
indexMin=min(find(fb >= p.fMin_Hz));
flag=1;
end
if isfield(p, 'fMax_Hz')
indexMax=max(find(fb <= p.fMax_Hz));
flag=1;
end
if flag
fb=fb(indexMin:indexMax);
fr=fr(indexMin:indexMax);
end

if p.NMax > 0
if p.logscale
% ************************************
% Sample number reduction for logarithmic
% frequency scale
% ************************************
assert(isfield(p, 'fMin_Hz'), 'need fMin_Hz in logscale mode');
assert(p.fMin_Hz > 0, 'need fMin_Hz > 0 in logscale mode');
if ~isfield(p, 'fMax_Hz')
p.fMax_Hz = p.rate_Hz / 2;
end

% averaged output arrays
fbOut=zeros(1, p.NMax-1);
frOut=zeros(1, p.NMax-1);

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

f1=s(1);
nextStartIndex=max(find(fb < f1));
if isempty(nextStartIndex)
nextStartIndex=1;
end

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

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

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

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

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

% ************************************
% Post-processing:
% Translate frequency axis to RF
% ************************************
fb = fb + p.originMapsTo_Hz;

% ************************************
% convert to requested units
% ************************************
switch(p.freqUnit)
case 'Hz'
case 'kHz'
fb = fb / 1e3;
case 'MHz'
fb = fb / 1e6;
case 'GHz'
fb = fb / 1e9;
otherwise
error('unsupported frequency unit');
end

% ************************************
% Plot (if requested)
% ************************************
if p.fig > 0
% *************************************************************
% title
% *************************************************************
if isfield(p, 'title')
t=['"', p.title, '";'];
else
t='';
end
switch(p.filter)
case 'gaussian'
t=[t, ' filter: Gaussian '];
case 'brickwall'
t=[t, ' filter: ideal bandpass '];
case 'none'
t=[t, ' filter: none '];
otherwise
assert(0)
end

if ~strcmp(p.filter, 'none')
t=[t, '(', mat2str(p.RBW_window_Hz), ' Hz)'];
end

if isfield(p, 'pNom_dBm')
% *************************************************************
% show dB relative to a given absolute power in dBm
% *************************************************************
fr=fr-p.pNom_dBm;
yUnit='dB';
else
yUnit='dBm';
end

% *************************************************************
% plot
% *************************************************************
figure(p.fig);
if p.logscale
handle = semilogx(fb, fr, p.plotStyle);
else
handle = plot(fb, fr, p.plotStyle);
end
hold on; % after plot, otherwise prevents logscale

if isfield(p, 'lineWidth')
set(handle, 'lineWidth', p.lineWidth);
end

% *************************************************************
% axis labels
% *************************************************************
xlabel(p.freqUnit);
if singleBinMode
ylabel([yUnit, ' (continuous wave)']);
else
ylabel([yUnit, ' in ', mat2str(p.RBW_power_Hz), ' Hz integration BW']);
end
title(t);

% *************************************************************
% adapt y range, if requested
% *************************************************************
y=ylim();
y1=y(1); y2=y(2);

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

function [winLeft, winRight] = spectrumAnalyzerGetWindow(filter, singleBinMode, RBW_window_Hz, binspacing_Hz)

switch(filter)
case 'gaussian'

% construct Gaussian filter
% -60 / -3 dB shape factor 4.472
nRBW=6;
nOneSide=ceil(RBW_window_Hz/binspacing_Hz*nRBW);

filterBase=linspace(0, nRBW, nOneSide);
winLeft=exp(-(filterBase*0.831) .^2);
winRight=fliplr(winLeft(2:end)); % omit 0 Hz bin

case 'brickwall'
nRBW=1;
n=ceil(RBW_window_Hz/binspacing_Hz*nRBW);
n1 = floor(n/2);
n2 = n - n1;

winLeft=ones(1, n1);
winRight=ones(1, n2);

case 'none'
winLeft=[];
winRight=[];

otherwise
error('unknown RBW filter type');
end

% the window is not supposed to shift the spectrum.
% Therefore, the first bin is always in winLeft (0 Hz):
if size(winLeft, 2)==1 && isempty(winRight)
% there is no use to convolve with one-sample window
% it's always unity
winLeft=[];
winRight=[];
tmpwin=[];
end

if ~isempty(winLeft)
% (note: it is not possible that winRight is empty, while winLeft is not)
if singleBinMode
% normalize to unity at 0 Hz
s=winLeft(1);
else
% normalize to unity area under the filter
s=sum(winLeft)+sum(winRight);
end
winLeft=winLeft / s;
winRight=winRight / s;
end
end

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

function r = flattenCellArray(arr, r)
ix=1;
ixMax = numel(arr);
while ix <= ixMax
e = arr{ix};

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

## Peak-to-average ratio analysis

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

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

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

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

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

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

assert(peakProbability > 0);
assert(peakProbability <= 1);

% *************************************************************
% determine RMS average and normalize signal to unity
% *************************************************************
nSamples = numel(signal);
sAverage = sqrt(signal * signal' / nSamples);
signal = signal / sAverage;

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

% *************************************************************
% determine number of samples with magnitude below the "peak"
% level, based on the given peak probability
% for example: probability = 0.5 => nBelowPeakLevel = nSamples/2
% typically, this will be a floating point number
% *************************************************************
nBelowPeakLevel = peakProbability * nSamples;

% *************************************************************
% interpolate the magnitude between the two closest samples
% *************************************************************
sMagSorted = sort(sMag);
x = [0 1:nSamples nSamples+1];
y = [0 sMagSorted sMagSorted(end)];
magAtPeakLevel = interp1(x, y, nBelowPeakLevel, 'linear');

% *************************************************************
% Peak-to-average ratio
% signal is normalized, average is 1
% the ratio relates to the sample magnitude
% "power" is square of magnitude => use 20*log10 for dB conversion.
% *************************************************************
PAR = magAtPeakLevel;
PAR_dB = 20*log10(PAR + eps);

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

## Phase continuity of NCO

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

test = 1;                              %see comments below
n = 10000;                             %number of test samples

switch test
case 1, z = linspace(.02,-.01,n);     %gentle zero crossing
case 2, z = linspace(.0005,-.0005,n); %excess dc at zero crossing
case 3, z = randsrc(1,n,[.25 -.25]);  %random sweep, creates messy signal
case 4, z = randsrc(1,n,[-.25:.001,.25]);  %random sweep, creates messy signal
case 5, z = ones(1,n/4)*.01; z = [z -z z -z]; %discontinuity at zero crossing
end

%%%%%%%%%%%%%%%%%%%%%%%%% Finite NCO model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 2^15 - 1;                          %max amplitude
M = 2^12;                              %NCO accumulator phase resolution
inc = round(M*z);                      %NCO accumulator phase increment
k = 2^12;                              %lut phase resolution (lut size),
lsb = log2(M) - log2(k);               %LSBs discarded when addressing
lut = round(A*exp(j*2*pi*(0:k-1)/k));  %lut, one cycle cos/sin data

ptr = 0;
addr = 0;
for i = 1:n
y(i) = lut(addr+1);                %add 1 for matlab LUT
ptr = mod(ptr + inc(i), M);        %phase accumulator(ramp)
addr = round(ptr/2^lsb);           %discard LSBs
addr(addr >= k) = addr - k;        %check address overflow
end

%display output
plot(real(y),'-');hold
plot(imag(y),'r-');``````

## How to speed up the sinc series

December 7, 2011 Coded in Matlab
``````function Demo

close all

T = 1; %Sampling period.
B = 0.7; %Two-sided bandwidth.

%Plot error spectrum for P = 10, 20, and 30.

PlotErrorSpectrum(T,B,10);

figure

PlotErrorSpectrum(T,B,20);

figure

PlotErrorSpectrum(T,B,30);

function PlotErrorSpectrum(T,B,P)

% Frequency grid specification in figure.
f0 = 0;
M = (2*P+1)*10;
Df = B/(2*M);

f = f0+Df*(0:M-1);

ESinc = zeros(size(f));
Eg = zeros(size(f));

%Look for the maximum error among different shifts u for each frequency in f.
for u = -0.5:1/((2*P+1)*5):0.5

HRef = exp(1i*2*pi*f*u); %Ideal spectrum.
HSinc = chirpzDFT(sinc((-P:P)+u/T),-P*T,T,f0,Df,M); %Sinc spectrum.
Hg = chirpzDFT(gPulseCoefs(u,T,B,P),-P*T,T,f0,Df,M); %g-pulse spectrum.

ESinc = max([ESinc;abs(HRef-HSinc)]); %Compute current max. error for sinc.
Eg = max([Eg;abs(HRef-Hg)]); %Compute current max. error for g.

end

plot(f,20*log10([ESinc;Eg].'))

xlabel('Freq. (f*T)')
ylabel('Maximum error (dB)')

legend('sinc pulse','g pulse','Location','NorthOutside')

pr = get(gca,'YTick');

text(0,2*pr(end)-pr(end-1),['B*T = ',num2str(B*T),', P = ',num2str(P)])

title('Error spectrum (maximum for all possible shifts u)')

grid on

function c = gPulseCoefs(u,T,B,P)

t = (-P:P)*T+u;
c = sinc(t/T).*real(sinc((1-B*T)*sqrt((t/T).^2-P^2)))/real(sinc(1i*(1-B*T)*P));

function X = chirpzDFT(x,to,Dt,fo,Df,M)

%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
%               N
%       X(k) = sum  x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
%              n=1
%
%Input parameters:
%
%  x: Signal samples.
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
%  M: Desired number of output samples.
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%

x = x(:).';

N = numel(x);
P = 2^ceil(log2(M+N-1));

%The next three lines do not depend on the vector x, and so they can be pre-computed if
%the time and frequency grids do not change, (i.e, for computing the transform of
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and
%three if the grids change.

a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse.

%Weigh x using a and perform FFT convolution with phi.
X = ifft( fft(x.*a,P) .* phi );

%Truncate the convolution tails and weigh using b.
X = X(N:M+N-1) .* b;``````

## I/Q and its conjugates

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

fc = .019;  % frequency relative to Fs of 1(-.5 ~ +.5)
phase = 0;  % degrees

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shift = floor(phase*1/(fc*360));
x1 = exp(j*2*pi*(0+shift:2^16-1+shift)*fc);   % I,Q
x2 = complex(real(x1),-imag(x1));             % I,-Q
x3 = complex(-real(x1),imag(x1));             %-I,Q
x4 = complex(imag(x1),real(x1));              % Q,I
x5 = complex(imag(x1),-real(x1));             % Q,-I
x6 = complex(-imag(x1),real(x1));             %-Q,I

f1 = fftshift(fft(x1,2^16));
f2 = fftshift(fft(x2,2^16));
f3 = fftshift(fft(x3,2^16));
f4 = fftshift(fft(x4,2^16));
f5 = fftshift(fft(x5,2^16));
f6 = fftshift(fft(x6,2^16));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = linspace(-.5,.5,2^16);
figure(1);
subplot(4,1,1);hold;
plot(f,20*log10(abs(f1)),'.-');
plot(f,20*log10(abs(f2)),'r-');
plot(f,20*log10(abs(f3)),'g--');
legend('I,Q','I,-Q','-I,Q')
ylabel('dB')

subplot(4,1,2);hold
plot(f,angle(f1),'.-');
plot(f,angle(f2),'r--');
plot(f,angle(f3),'g--');
legend('I,Q','I,-Q','-I,Q')
ylabel('rad')

subplot(4,1,3);hold;
plot(f,20*log10(abs(f4)));
plot(f,20*log10(abs(f5)),'r-');
plot(f,20*log10(abs(f6)),'g--');
legend('Q,I','Q,-I','-Q,I')
ylabel('dB')

subplot(4,1,4);hold
plot(f,angle(f4));
plot(f,angle(f5),'r--');
plot(f,angle(f6),'g--');
legend('Q,I','Q,-I','-Q,I')
ylabel('rad')

z = 1:ceil(1/abs(fc));
figure(2)
subplot(6,1,1);hold
plot(real(x1(z)),'.-');
plot(imag(x1(z)),'r.-');
legend('I','Q');

subplot(6,1,2);hold
plot(real(x2(z)),'.-');
plot(imag(x2(z)),'r.-');
legend('I','-Q');

subplot(6,1,3);hold
plot(real(x3(z)),'.-');
plot(imag(x3(z)),'r.-');
legend('-I','Q');

subplot(6,1,4);hold
plot(real(x4(z)),'.-');
plot(imag(x4(z)),'r.-');
legend('Q','I');

subplot(6,1,5);hold
plot(real(x5(z)),'.-');
plot(imag(x5(z)),'r.-');
legend('Q','-I');

subplot(6,1,6);hold
plot(real(x6(z)),'.-');
plot(imag(x6(z)),'r.-');
legend('-Q','I');``````

## Tremolo audio effect (amplitude modulation)

November 29, 20113 comments Coded in C
``````/************************Tremolo.c*************************/

#include "Tremolo.h"

static double dep;
static short counter_limit;
static short control;
static short mod;
static double offset;

void Tremolo_init(short effect_rate,double depth) {
dep     = depth;
control = 1;
mod     = 0;
counter_limit = effect_rate;
offset  = 1 - dep;
}

double Tremolo_process(double xin) {
double yout;
double m;

m = (double)mod*dep/counter_limit;
yout = (m + offset)*xin;
return yout;
}

void Tremolo_sweep(void) {

mod += control;

if (mod > counter_limit) {
control = -1;
}
else if(!mod) {
control = 1;
}
}

/********************Tremolo.h****************************/
#ifndef __TREMOLO_H__
#define __TREMOLO_H__

extern void Tremolo_init(short effect_rate,double depth);
extern double Tremolo_process(double xin);
extern void Tremolo_sweep(void);

#endif

/*Usage example*/
#include "Tremolo.h"

void main(void) {
double xin;
double yout;
Tremolo_init(4000,1);

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

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

## Auto Wah and Envelope follower audio effect

November 29, 2011 Coded in C
``````#include "bp_iir.h"
#include "AutoWah.h"

static short center_freq;
static short samp_freq;
static short counter;
static short counter_limit;
static short control;
static short max_freq;
static short min_freq;
static short f_step;
static struct bp_filter H;
static double a[3];
static double b[3];
double x[3],y[3];

void AutoWah_init(short effect_rate,short sampling,short maxf,short minf,short Q,double gainfactor,short freq_step) {
double C;

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

//User Parametters
counter_limit = effect_rate;

//Convert frequencies to index ranges
min_freq = 0;
max_freq = (maxf - minf)/freq_step;

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

//Lowpass filter parameters
/*
C = 1/TAN(PI*Fc/Fs)
*/
C = 1018.59;

b[0] = 1/(1+2*0.7071*C+pow(C,2));
b[1] = 2*b[0];
b[2] = b[0];

a[1] = 2*b[0]*(1-pow(C,2));
a[2] = b[0]*(1-2*0.7071*C+pow(C,2));
}

double AutoWah_process(double xin) {
double yout;

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

return yout;
}

void AutoWah_sweep(double xin) {
unsigned int filter_index;
double yout;
double detect;

/*The input is 16 bit unsigned so it
has to be centered to 0*/
detect = (xin - 32767.0);
x[0] = x[1];
x[1] = x[2];
/*Here the input to the filter
is half rectified*/
x[2] = (detect > 0)?detect:0;

y[0] = y[1];
y[1] = y[2];

y[2] =    b[0]*x[2];
y[2] +=   b[1]*x[1];
y[2] +=   b[2]*x[0];
y[2] -=   a[1]*y[1];
y[2] -=   a[2]*y[0];

if (!--counter) {
/*The output of the LPF (y[2]) is scaled by 0.1
in order to be used as a LFO to control the band pass filter
*/
filter_index = (double)min_freq + y[2]*0.1;

/*The scaling value is determined as a value
that would keep the filter index within the
range of max_freq
*/
if(filter_index > max_freq) {
filter_index = max_freq;
}
if(filter_index < 0) {
filter_index = 0;
}

bp_iir_setup(&H,filter_index);

counter = counter_limit;
}
}``````

## Computing FFT Twiddle Factors

November 28, 201110 comments Coded in Matlab
``````% Filename:  FFT_Twiddles_Find_DSPrelated.m

%  Computes 'Decimation in Frequency' or 'Decimation
%  in Time' Butterfly twiddle factors, for radix-2 FFTs
%  with in-order input indices and scrambled output indices.
%
%  To use, do two things: (1) define FFT size 'N'; and
%  (2) define the desired 'Structure', near line 17-18,
%  as 'Dec_in_Time' or 'Dec_in_Freq'.
%
%  Author: Richard Lyons, November, 2011
%******************************************

clear, clc

% Define input parameters
N = 8; % FFT size (Must be an integer power of 2)
Structure = 'Dec_in_Time';  % Choose Dec-in-time butterflies
%Structure = 'Dec_in_Freq';  % Choose Dec-in-frequency butterflies

% Start of processing
Num_Stages = log2(N); % Number of stages
StageStart = 1; % First stage to compute
StageStop = Num_Stages; % Last stage to compute
ButterStart = 1; %First butterfly to compute
ButterStop = N/2; %Last butterfly to compute
Pointer = 0; %Init 'results' row pointer
for Stage_Num = StageStart:StageStop
if Structure == 'Dec_in_Time'
for Butter_Num = ButterStart:ButterStop
Twid = floor((2^Stage_Num*(Butter_Num-1))/N);
% Compute bit reversal of Twid
Twid_Bit_Rev = 0;
for I = Num_Stages-2:-1:0
if Twid>=2^I
Twid_Bit_Rev = Twid_Bit_Rev + 2^(Num_Stages-I-2);
Twid = Twid -2^I;
else, end
end %End bit reversal 'I' loop
A1 = Twid_Bit_Rev; %Angle A1
A2 = Twid_Bit_Rev + N/2; %Angle A2
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Butter_Num,A1,A2];
end
else
for Twiddle_Num = 1:N/2^Stage_Num
Twid = (2^Stage_Num*(Twiddle_Num-1))/2; %Compute integer
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Twiddle_Num,Twid];
end
end % End 'if'
end % End Stage_Num loop

Results(:,1:3), disp(' Stage#  Twid#   A'), disp(' ')``````

## Computing Acceptable Bandpass Sample Rates

November 25, 2011 Coded in Matlab
``````% Filename: Bandpass_Sample_Rate_Calc.m
%
%  Calculates acceptable Bandpass Sampling rate ranges
%  based on an analog (continuous) signal's bandwidth
%  and center freq.
%
%  Merely define bandwidth "Bw" and center frequency "Fc", in
%  Hz, near line 22, for the analog bandpass signal and run the
%  program.  [Example: Bw = 5, Fc = 20.]  Acceptable Fs sample
%  rate ranges are shown in Figure 1 and displayed in Command window.
%
%  If the User defines a value for the BP sample rate Fs, near
%  near line 28, then Figure 2 will show the locations of the
%  positive and negative-freq spectral components after
%  bandpass sampling.
%
%   Richard Lyons [November 2011]
%******************************************

clear, clc

Bw = 5; % Analog signal bandwidth in Hz
Fc = 20; % Analog signal center freq in Hz

% ##############################################
% Define an Fs sample rate value below

Fs = 11; % Selected Fs sample rate in Hz
% ##############################################

disp(' '), disp(['Analog Center Freq = ',num2str(Fc),' Hz.'])
disp(['Analog Bandwidth = ',num2str(Bw),' Hz.']), disp(' ')

% *****************************************************
%  Compute % display the acceptable ranges of BP sample rate
% *****************************************************
disp('----------------------------------')
disp('Acceptable Fs ranges in Hz:')
No_aliasing = 0;  % Init a warning flag
M = 1; % Initialize a counter

while (2*Fc + Bw)/(M+1) >= 2*Bw
FsMin = (2*Fc + Bw)/(M+1);
FsMax = (2*Fc - Bw)/M;
Fs_ranges(M,1) = FsMin;
Fs_ranges(M,2) = FsMax;
M = M + 1;
disp([num2str(FsMin),' -to- ',num2str(FsMax)])
end
disp('----------------------------------')

% *****************************************************
%  Plot the acceptable ranges of BP sample rate
% *****************************************************
figure(1), clf
title('Acceptable Ranges of Bandpass Sampling Rate')
xlabel('Freq (Hz)')
ylabel('This axis has no meaning')

for K = 1:M-1
hold on
plot([Fs_ranges(K,1),Fs_ranges(K,1)],[0.5,1.2],':g');
plot([Fs_ranges(K,1),Fs_ranges(K,2)],[1,1],'-r','linewidth',4);
axis([0.8*(2*Bw),1.1*max(Fs_ranges(1,2)), 0.8, 1.55])
end

plot([2*Bw,2*Bw],[0.5,1.2],'-b','linewidth',2);
text(2*Bw,1.45,'Bold red lines show acceptable Fs ranges')
text(2*Bw,1.35,'Blue line = Twice analog signal Bandwidth (2 x Bw)')
text(2*Bw,1.25,'(You can zoom in, if you wish.)')
hold off, grid on, zoom on

% **************************************************************
%  If Fs has been defined, plot spectrum of the sampled signal.
% **************************************************************
%
% Check if "Fs" has been defined
disp(' ')
if isempty(strmatch('Fs',who,'exact')) == 1;
disp('Fs sampling rate has NOT been defined.')
% Fs does NOT exist, do nothing further.
else
% Fs is defined, plot the spectrum of the sampled signal.
disp(['Fs defined as ',num2str(Fs),' Hz'])

% To determine intermediate frequency (IF), check integer
%      part of "2Fc/Fs" for odd or even
Temp = floor(2*Fc/Fs);
if (Temp == 2*floor(Temp/2))
disp(' '), disp('Pos-freq sampled spectra is not inverted.')
Freq_IF = Fc -Fs*floor(Fc/Fs); % Computed IF frequency
else
disp(' '), disp('Pos-freq sampled spectra is inverted.')
Freq_IF = Fs*(1 + floor(Fc/Fs)) -Fc; % Computed IF frequency
end
disp(' '), disp(['Center of pos-freq sampled spectra = ',num2str(Freq_IF),' Hz.'])

% Prepare to plot sampled spectral range in Figure 2
IF_MinFreq = Freq_IF-Bw/2;
IF_MaxFreq = Freq_IF+Bw/2;

figure(2), clf
hold on
plot([IF_MinFreq,IF_MaxFreq],[0.95, 1],'-r','linewidth',4);
plot([Fs-IF_MaxFreq,Fs-IF_MinFreq],[1, 0.95],'-r','linewidth',4);
plot([Fs,Fs],[0.5, 1.02],'-b','linewidth',2);
plot([Fs/2,Fs/2],[0.5, 1.02],':b','linewidth',2);
plot([IF_MinFreq,IF_MinFreq],[0.5, 0.95],':r');
plot([IF_MaxFreq,IF_MaxFreq],[0.5, 1],':r');
plot([Fs-IF_MaxFreq,Fs-IF_MaxFreq],[0.5, 1],':r');
plot([Fs-IF_MinFreq,Fs-IF_MinFreq],[0.5, 0.95],':r');
text(0.9*Fs,1.03,['Fs = ',num2str(Fs),' Hz'])
text(0.8*Fs/2, 1.03,['Fs/2 = ',num2str(Fs/2),' Hz'])
text(Fs/10,1.07,'(You can zoom in, if you wish.)')
axis([0,1.2*Fs, 0.8, 1.1])

hold off
title('Red lines show spectral range of sampled signal')
xlabel('Freq (Hz)')
ylabel('This axis has no meaning')
grid on, zoom on

% ################################################################
% Check if Fs is NOT in an acceptable freq range (aliasing occurs)
Aliasing_Flag = 1; % Initialize a flag
for K = 1:M-1 % Check each individual acceptable Fs range
if Fs_ranges(K,1)<=Fs & Fs<=Fs_ranges(K,2)% & Fs<=Fs_ranges(K,2)
% No aliasing will occur
Aliasing_Flag = Aliasing_Flag -1;
else, end
end % End K loop
if Aliasing_Flag == 1;  % Aliasing will occur
text(Fs/10, 0.91, 'WARNING! WARNING!')
text(Fs/10, 0.89, 'Aliasing will occur!')
else, end
zoom on
% ################################################################
end``````