## Peak-to-average ratio analysis

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

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

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

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

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

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

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

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

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

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

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

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

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

## 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),
lut = round(A*exp(j*2*pi*(0:k-1)/k));  %lut, one cycle cos/sin data

ptr = 0;
for i = 1:n
ptr = mod(ptr + inc(i), M);        %phase accumulator(ramp)
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')

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

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

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

## Computing Optimum Two-Stage Decimation Factors

November 21, 2011 Coded in Matlab
function [D1,D2] = Dec_OptimTwoStage_Factors(D_Total,Passband_width,Fstop)
%
%    When breaking a single large decimation factor (D_Total) into
%    two stages of decimation (to minimize the orders of the
%    necessary lowpass digital filtering), an optimum first
%    stage of decimation factor (D1) can be found.  The second
%    stage's decimation factor is then D_Total/D1.
%
%    Inputs:
%         D_Total = original single-stage decimation factor.
%         Passband_width = desired passband width of original
%                          single-stage lowpass filter (Hz).
%         Fstop = desired beginning of stopband freq of original
%                          single-stage lowpass filter (Hz).
%
%    Outputs:
%         D1 = optimum first-stage decimation factor.
%         D2 = optimum second-stage decimation factor.
%
%    Example: Assume you want to decimate a sequence by D_Total=100,
%             your original single-stage lowpass filter's passband
%             width and stopband frequency are 1800 and 2200 Hz
%             respectively.  We use:
%
%              [D1,D2] = Dec_OptimTwoStage_Factors(100, 1800, 2200)
%
%             giving us the desired D1 = 25, and D2 = 4.
%             (That is, D1*D2 = 25*4 = 100 = D_Total.)
%
%    Author: Richard Lyons, November, 2011

% Start of processing
DeltaF = (Fstop -Passband_width)/Fstop;

Radical = (D_Total*DeltaF)/(2 -DeltaF); % Terms under sqrt sign.
Numer = 2*D_Total*(1 -sqrt(Radical));  % Numerator.
Denom = 2 -DeltaF*(D_Total + 1);  % Denominator.
D1_estimated = Numer/Denom;  %Optimum D1 factor, but not
%                  an integer.

% Find all factors of 'D_Total'
Factors_of_D_Total = Find_Factors_of_D_Total(D_Total);

% Find the one factor of 'D_Total' that's closest to 'D1_estimated'
Temp = abs(Factors_of_D_Total -D1_estimated);
Index_of_Min = find(Temp == min(Temp)); % Index of optim. D1
D1 = Factors_of_D_Total(Index_of_Min); % Optimum first
%    decimation factor
D2 = D_Total/D1;  % Second decimation factor.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allfactors] = Find_Factors_of_D_Total(X)
% Finds all integer factors of intger X.
% Filename was originally 'listfactors.m', written
[b,m,n] = unique(factor(X));
%b is all prime factors
occurences = [m(1) diff(m)];
current_factors = [b(1).^(0:occurences(1))]';
for index = 2:length(b)
currentprime = b(index).^(0:occurences(index));
current_factors = current_factors*currentprime;
current_factors = current_factors(:);
end
allfactors = sort(current_factors);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## FFT for arbitrary time and frequency grids (Chirp-z Transform)

November 17, 2011 Coded in Matlab
function Demo

while 1

%Frequency grid with arbitrary start value, spacing, and number of elements. ...
%Mathematically, the grid is
%
%   fo+Df*m, m=0, 1,...,M-1.
%

fo = 10*rand(1)-0.5; %Start value.
Df = rand(1)+0.01;   %Frequency spacing.
M = round(rand(1)*1000); %Number of frequencies.

%Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
%the grid is
%
%   to+Dt*n, n=0, 1,...,N-1.
%

to = 10*rand(1)-0.5; %Start value; (instant of first sample).
Dt = rand(1)+0.01;   %Time spacing.
N = round(rand(1)*1000); %Number of samples.

%Random vector of samples.
x = randn(1,N)+1i*randn(1,N);

%x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N.

%We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.

%Compute DFT using the direct and chirp-z methods.

tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;

disp(['Timing direct method: ', num2str(tmD), ' sec.'])
disp(['Timing chirp z: ',num2str(tmF),' sec.'])
disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
disp(' ')
input('(Press a key for another trial or Ctrl-C) ')
disp(' ')
end

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;

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

%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.

x = x(:).';

N = length(x);

X = zeros(1,M);

for k = 1:M
for n = 1:N
X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
end
end

## NCO bit true model in a single equation

October 23, 2011 Coded in Matlab
%NCO bit true model
clear all; close all;

%example NCO parameters
n = 20000;                     %number of test samples
A = 2^15 - 1;                  %max amplitude, amplitude resolution 16 bit signed
M = 2^32;                      %NCO accumulator phase resolution
inc = round(M*.001);           %NCO accumulator phase increment
k = 2^12;                      %lut phase resolution (lut size), may equal M

%single equation for NCO bit true model
y1 = round(A*sin(2*pi*round(k*mod((0:n-1)*inc/M,M))/k));

%%%%%%%%%%%%%%%%%%%%%%%%% equivalent code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lut = round(A*sin(2*pi*(0:k-1)/k));      %lut, one cycle sine data

ptr = 0;
for i = 1:n
ptr = mod(ptr + inc, M);      %phase accumulator(ramp)
end

%compare
plot(y1);hold;
plot(y2,'r--')
plot(y1 - y2,'g--')

## Complex filter with three multipliers per tap

October 22, 2011 Coded in Matlab
%random input
I_in = randn(1,1024);
Q_in = randn(1,1024);

%random test filter
h = randn(1,20)+j*randn(1,20);

%split up into 3 subfilters
h2 = real(h);
h3 = real(h)-imag(h);         %extra subtractor

%apply subfilter
ch1 = filter(h1,1,I_in);
ch2 = filter(h2,1,Q_in-I_in); %extra subtractor
ch3 = filter(h3,1,Q_in);

%combine outputs
yr = round((ch3-ch2)*2^15);
yi = round((ch1+ch2)*2^15);
IQ_out = complex(yr,yi);

%direct computation
ref = round(2^15*filter(h,1,complex(I_in,Q_in)));
plot(IQ_out - ref);