Code Snippets Submitted by Rick Lyons

## Computing the Nth Roots of a Number

```
function [nth_Roots] = roots_nth(Num, n, Plot)
% ROOTS_NTH computes the nth roots of Num.
%
% Inputs:
% Num is the number(s) for which the nth roots will be computed.
% Num can be real-valued or complex-valued.
% If Num has more than one element, then Num must be a row vector.
% n must be a positive integer.
% Plot is a string equal to 'Y' if a complex plane plot of
% the n complex roots of each element in Num is desired.
% Principle root is plotted as a blue square. Remaining roots
% are red dots. If input Plot string does not exist, then no
% plot is generated.
%
% Output:
% nth_Roots is a column vector (with n rows) giving the nth roots of Num.
%
% Example:
% clear, clc % clear all variables
% Num = [2+j*11, 3*exp(j*pi)]; % Two numbers
% n = 3; % find the 3rd roots
% [nth_Roots] = roots_nth(Num, n, 'Y')
%
% Results are:
% nth_Roots =
%
% 2.0000 + 1.0000i 0.7211 + 1.2490i
% -1.8660 + 1.2321i -1.4422 + 0.0000i
% -0.1340 - 2.2321i 0.7211 - 1.2490i
%
% [R. Lyons, Jan. 2013]
% Input argument check
if (nargin<2),
disp('Not enough input arguments!!')
return
end;
if (nargin==2)
Plot == 'N'
end
Num_of_Columns = max(size(Num)); % How many elements are in Num?
for Column_Loop = 1:Num_of_Columns;
Mag(Column_Loop) = abs(Num(Column_Loop));
Angle_Degrees(Column_Loop) = angle(Num(Column_Loop))*180/pi;
for k = 1:n
Root_Mag = Mag(Column_Loop)^(1/n);
Root_Angle_Deg(k) = Angle_Degrees(Column_Loop)/n + (k-1)*360/n;
nth_Roots(k,Column_Loop) = Root_Mag*exp(j*Root_Angle_Deg(k)*pi/180);
end
% Plot roots, on complex plane(s), if necessary
if Plot == 'Y'
figure
% Plot unit circle
Num_Pts = 200; % Number of circle points
Index = 0 : Num_Pts-1;
Alpha = 1/Num_Pts;
Points = Root_Mag*exp(j*2*pi*Alpha*Index);
plot(Points, ':k','markersize',5)
grid on
hold on
% Plot principle root (square blue dot) and list results in plot
axis([-1.1*Root_Mag, 1.1*Root_Mag, -1.1*Root_Mag, 1.1*Root_Mag]);
plot(Root_Mag*cos(Root_Angle_Deg(1)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(1)*pi/180), 'bs');
text(-Root_Mag/2, 3*Root_Mag/4, ['Magnitudes = ',num2str(Root_Mag)]);
text(-Root_Mag/2, 3*Root_Mag/4-Root_Mag/10, ...
['Principle Angle (deg.) = ',num2str(Root_Angle_Deg(1))]);
% Plot remaining roots as red dots
for k = 2:n
plot(Root_Mag*cos(Root_Angle_Deg(k)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(k)*pi/180), 'ro');
text(-Root_Mag/2, 3*Root_Mag/4-k*Root_Mag/10, ...
[' Angle ',num2str(k),' (deg.) = ', ...
num2str(Root_Angle_Deg(k))]);
end % End k loop
xlabel('Real'); ylabel('Imag.');
hold off
end % End 'if Plot == 'Y''
end % End Column_Loop
```

## Computing FFT Twiddle Factors

```
% 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(' ')
```

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

```
function [Windowed_Spec] = Wind_Flattop(Spec)
% Given an input spectral sequence 'Spec', that is the
% FFT of some time sequence 'x', Wind_Flattop(Spec)
% returns a spectral sequence that is equivalent
% to the FFT of a flat-top windowed version of time
% sequence 'x'. The peak magnitude values of output
% sequence 'Windowed_Spec' can be used to accurately
% estimate the peak amplitudes of sinusoidal components
% in time sequence 'x'.
% Input: 'Spec' (a sequence of complex FFT sample values)
%
% Output: 'Windowed_Spec' (a sequence of complex flat-top
% windowed FFT sample values)
%
% Based on Lyons': "Reducing FFT Scalloping Loss Errors
% Without Multiplication", IEEE Signal Processing Magazine,
% DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
% Richard Lyons [December, 2011]
N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];
% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);
Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);
% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);
Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);
% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function
```

## Computing Acceptable Bandpass Sample Rates

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

## Beat Notes

```
% Filename: Beat_Frequency.m
%
% [Richard Lyons, Feb. 2013]
clear, clc
Fs = 8192; % Sample rate of dig. samples
N = 8192; % Number of time samples
n = 0:N-1;
Wave_1 = sin(2*pi*210*n/Fs); % First tone, 210 Hz
Wave_2 = sin(2*pi*200*n/Fs); % Second tone, 200 Hz
% Plot the two tones
figure(1)
plot(n/Fs, Wave_1, '-b')
ylabel('200 Hz'); xlabel('Time (sec.)');
hold on
plot(n/Fs, Wave_2, '-r')
axis([0, 0.05, -1.2, 1.5]);
ylabel('Input tones'); xlabel('Time (sec.)');
title('red = 200 Hz tone, blue = 210 Hz tone');
grid on, zoom on
hold off
Product = Wave_1.*Wave_2;
Sum = Wave_1 + Wave_2;
% Plot the tones' product and sum
figure(2)
subplot(2,1,1)
plot(n/Fs, Product, '-b'),
ylabel('Product'); xlabel('Time (sec.)');
grid on, zoom on
hold on
Red_Curve = 0.5*cos(2*pi*10*n/Fs) + 0.5; % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -1.25, 1.5]);
hold off
grid on, zoom on
subplot(2,1,2)
plot(n/Fs, Sum, '-b')
hold on
Red_Curve = 2*cos(2*pi*5*n/Fs); % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -2.4, 3]);
hold off
ylabel('Sum'); xlabel('Time (sec.)');
grid on, zoom on
% Play all the signals
sound(Wave_1, Fs)
pause(1.2)
sound(Wave_2, Fs)
pause(1.2)
sound(Product, Fs)
pause(1.2)
sound(Sum, Fs)
% Spec analysis of the "Sum" signal
Spec = fft(Sum);
Spec_Mag = abs(Spec);
Freq = n*Fs/N; % Freq vector in Hz
figure (3) % Plot positive-freq spec amg
plot(Freq(1:N/16), Spec_Mag(1:N/16))
title('Spec Mag of "Sum" signal')
ylabel('Mag.'), xlabel('Hz')
grid on, zoom on
```

## Testing the Flat-Top Windowing Function

```
% Code for testing the 'Wind_Flattop(Spec)' function in
% reducing 'scalloping loss' errors in time signal amplitude
% estimation.
%
% Generates a time-domain sinusoid, computes its FFT, and passes
% that FFT sequence to the 'Wind_Flattop(Spec)' function.
%
% The maximum output sample of the 'Wind_Flattop(Spec)' function
% is used to estimate the peak amplitude of the original
% sinusoidal time-domain test signal.
%
% The user controls the values for the test sinusoid's
% 'Test_Freq' and peak amplitude 'Peak_Amp' near lines 23 & 24.
%
% Richard Lyons [December, 2011]
clear all, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define test parameters
Test_Freq = 7.22; % Test tone's freq. Must be less that N/2
Peak_Amp = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 64; % Number of time samples
Index = (0:N-1);
X = Peak_Amp*cos(2*pi*(Test_Freq)*Index/N + pi/3);
figure(1), clf
subplot(2,1,1)
plot(X,':ko', 'markersize', 4)
title('Original time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT the input 'X' sequence and call 'Wind_Flattop()' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Spec = fft(X);
[Windowed_Spec] = Wind_Flattop(Spec);
subplot(2,1,2), plot(abs(Spec),'ko', 'markersize', 3)
title('SpecMag of unwindowed time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results accuracy (error)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(['Test Freq = ',num2str(Test_Freq),...
', True Peak Amplitude = ',num2str(Peak_Amp)])
Mag_peak_unwindowed = max(abs(Spec));
Unwindowed_Amp_Estimate = 2*Mag_peak_unwindowed/N;
Unwindowed_Amp_Estimate_Error_in_dB = ...
20*log10(Unwindowed_Amp_Estimate/Peak_Amp);
disp(' ')
disp(['Unwindowed Peak Amplitude Estimate = ',...
num2str(Unwindowed_Amp_Estimate)])
disp(['Unwindowed Estimate Error in dB = ',...
num2str(Unwindowed_Amp_Estimate_Error_in_dB),' dB'])
M_peak_windowed = max(abs(Windowed_Spec));
Windowed_Amp_Estimated = 2*M_peak_windowed/N;
Windowed_Amp_Estimation_Error_in_dB = ...
20*log10(Windowed_Amp_Estimated/Peak_Amp);
disp(' ')
disp(['Windowed Peak Amplitude Estimate = ',...
num2str(Windowed_Amp_Estimated)])
disp(['Windowed Estimate Error in dB = ',...
num2str(Windowed_Amp_Estimation_Error_in_dB),' dB'])
```

## Computing Optimum Two-Stage Decimation Factors

```
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
% by Jeff Miller. Downloaded from Matlab Central.
[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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

## Computing CIC Filter Register Pruning Using Matlab [Updated]

```
% Filename: CIC_Word_Truncation.m
%
% Computes CIC decimation filter accumulator register
% truncation in each filter stage based on Hogenauer's
% 'accumulator register pruning' technique.
%
% Inputs:
% N = number of decimation CIC filter stages (filter order).
% R = CIC filter rate change factor (decimation factor).
% M = differential delay.
% Bin = number of bits in an input data word.
% Bout = number of bits in the filter's final output data word.
% Outputs:
% Stage number (ranges from 1 -to- 2*N+1).
% Bj = number of least significant bits that can be truncated
% at the input of a filter stage.
% Accumulator widths = number of a stage's necessary accumulator
% bits accounting for truncation.
% Richard Lyons Feb., 2012
clear, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define CIC filter parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N = 4; R = 25; M = 1; Bin = 16; Bout = 16; % Hogenauer paper, pp. 159
%N = 3; R = 32; M = 2; Bin = 8; Bout = 10; % Meyer Baese book, pp. 268
%N = 3; R = 16; M = 1; Bin = 16; Bout = 16; % Thorwartl's PDF file
N = 3; R = 8; M = 1; Bin = 12; Bout = 12; % Lyons' blog Figure 2 example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find h_sub_j and "F_sub_j" values for (N-1) cascaded integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
for j = N-1:-1:1
h_sub_j = [];
h_sub_j((R*M-1)*N + j -1 + 1) = 0;
for k = 0:(R*M-1)*N + j -1
for L = 0:floor(k/(R*M)) % Use uppercase "L" for loop variable
Change_to_Result = (-1)^L*nchoosek(N, L)...
*nchoosek(N-j+k-R*M*L,k-R*M*L);
h_sub_j(k+1) = h_sub_j(k+1) + Change_to_Result;
end % End "L" loop
end % End "k" loop
F_sub_j(j) = sqrt(sum(h_sub_j.^2));
end % End "j" loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for up to seven cascaded combs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j_for_many_combs = sqrt([2, 6, 20, 70, 252, 924, 3432]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for last integrator stage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(N) = F_sub_j_for_many_combs(N-1)*sqrt(R*M); % Last integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for N cascaded filter's comb stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2*N:-1:N+1
F_sub_j(j) = F_sub_j_for_many_combs(2*N -j + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for the final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(2*N+1) = 1; % Final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of minus log base 2 of "F_sub_j" values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minus_log2_of_F_sub_j = -log2(F_sub_j)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute total "Output_Truncation_Noise_Variance" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CIC_Filter_Gain = (R*M)^N;
Num_of_Bits_Growth = ceil(log2(CIC_Filter_Gain));
% The following is from Hogenauer's Eq. (11)
%Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin -1;
Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin;
Num_of_Output_Bits_Truncated = Num_Output_Bits_With_No_Truncation -Bout;
Output_Truncation_Noise_Variance = (2^Num_of_Output_Bits_Truncated)^2/12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute log base 2 of "Output_Truncation_Noise_Standard_Deviation" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output_Truncation_Noise_Standard_Deviation = ...
sqrt(Output_Truncation_Noise_Variance);
Log_base2_of_Output_Truncation_Noise_Standard_Deviation = ...
log2(Output_Truncation_Noise_Standard_Deviation);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of "half log base 2 of 6/N" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Half_Log_Base2_of_6_over_N = 0.5*log2(6/N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute desired "B_sub_j" vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_sub_j = floor(Minus_log2_of_F_sub_j ...
+ Log_base2_of_Output_Truncation_Noise_Standard_Deviation ...
+ Half_Log_Base2_of_6_over_N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '), disp(' ')
disp(['N = ',num2str(N),', R = ',num2str(R),', M = ',num2str(M),...
', Bin = ', num2str(Bin),', Bout = ',num2str(Bout)])
disp(['Num of Bits Growth Due To CIC Filter Gain = ', ...
num2str(Num_of_Bits_Growth)])
disp(['Num of Accumulator Bits With No Truncation = ', ...
num2str(Num_Output_Bits_With_No_Truncation)])
% disp(['Output Truncation Noise Variance = ', ...
% num2str(Output_Truncation_Noise_Variance)])
% disp(['Log Base2 of Output Truncation Noise Standard Deviation = ',...
% num2str(Log_base2_of_Output_Truncation_Noise_Standard_Deviation)])
% disp(['Half Log Base2 of 6/N = ', num2str(Half_Log_Base2_of_6_over_N)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and display "Results" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Stage = 1:2*N
Results(Stage,1) = Stage;
Results(Stage,2) = F_sub_j(Stage);
Results(Stage,3) = Minus_log2_of_F_sub_j(Stage);
Results(Stage,4) = B_sub_j(Stage);
Results(Stage,5) = Num_Output_Bits_With_No_Truncation -B_sub_j(Stage);
end
% Include final output stage truncation in "Results" matrix
Results(2*N+1,1) = 2*N+1; % Output stage number
Results(2*N+1,2) = 1;
Results(2*N+1,4) = Num_of_Output_Bits_Truncated;
Results(2*N+1,5) = Bout;
%Results % Display "Results" matrix in raw float-pt.form
% % Display "F_sub_j" values if you wish
% disp(' ')
% disp(' Stage Fj -log2(Fj) Bj Accum width')
% for Stage = 1:2*N+1
% disp([' ',sprintf('%2.2g',Results(Stage,1)),sprintf('\t'),sprintf('%12.3g',Results(Stage,2)),...
% sprintf('\t'),sprintf('%7.5g',Results(Stage,3)),sprintf('\t'),...
% sprintf('%7.5g',Results(Stage,4)),sprintf('\t'),sprintf('%7.5g',Results(Stage,5))])
% end
% Display Stage number, # of truncated input bits, & Accumulator word widths
disp(' ')
disp(' Stage(j) Bj Accum (adder) width')
for Stage = 1:2*N
disp([' ',sprintf('%2.0f',Results(Stage,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(Stage,4)),sprintf('\t'),...
sprintf('%7.5g',Results(Stage,5))])
end
disp([' ',sprintf('%2.0f',Results(2*N+1,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(2*N+1,4)),sprintf('\t'),...
sprintf('%7.5g',Results(2*N+1,5)),' (final truncation)'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

## Adding a Controlled Amount of Noise to a Noise-Free Signal

```
function [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
% SNR_Set(x, Desired_SNR_dB) returns the real-valued
% input 'Signal' contaminated with normally-distributed,
% zero-mean, random noise. The signal-to-noise ratio
% (SNR in dB) of the output 'Noisy_Signal' signal is
% controlled by the input 'Desired_SNR_dB' argument measured
% in dB.
% Example:
% Npts = 128; % Number of time samples
% n = 0:Npts-1; % Time-domain index
% Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% Desired_SNR_dB = 3; % Set SNR of output 'Noisy_Signal' to +3 dB
% [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
% Author: Richard Lyons [December 2011]
%******************************************
Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one
Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));
K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10); % Scale factor
New_Noise = sqrt(K)*Noise; % Change Noise level
%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))
Noisy_Signal = Signal + New_Noise;
'SNR_Set()' Function Test Code:
% Filename SNR_Set_test.m
%
% Tests the 'SNR_Set()" function. Adds a predefined
% amount of random noise to a noise-free signal such that
% the noisy signal has a desired signal-to-noise ratio (SNR).
%
% Author: Richard Lyons [December 2011]
clear, clc
% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5; % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal
Desired_SNR_dB = 3 % Set desired SNR in dB
[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal
% Plot original and 'noisy' signals
figure(1)
subplot(2,1,1)
plot(n, Signal, '-bo', 'markersize', 2)
title('Original Signal')
grid on, zoom on
subplot(2,1,2)
plot(n, Noisy_Signal, '-bo', 'markersize', 2)
title('Noisy Signal')
xlabel('Time-samples')
grid on, zoom on
% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
figure(2)
plot(Spec_Mag, '-bo', 'markersize', 2)
title('Spec Mag of Noisy Signal')
xlabel('Freq-samples'), ylabel('Linear')
grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)
```

## Beat Notes

```
% Filename: Beat_Frequency.m
%
% [Richard Lyons, Feb. 2013]
clear, clc
Fs = 8192; % Sample rate of dig. samples
N = 8192; % Number of time samples
n = 0:N-1;
Wave_1 = sin(2*pi*210*n/Fs); % First tone, 210 Hz
Wave_2 = sin(2*pi*200*n/Fs); % Second tone, 200 Hz
% Plot the two tones
figure(1)
plot(n/Fs, Wave_1, '-b')
ylabel('200 Hz'); xlabel('Time (sec.)');
hold on
plot(n/Fs, Wave_2, '-r')
axis([0, 0.05, -1.2, 1.5]);
ylabel('Input tones'); xlabel('Time (sec.)');
title('red = 200 Hz tone, blue = 210 Hz tone');
grid on, zoom on
hold off
Product = Wave_1.*Wave_2;
Sum = Wave_1 + Wave_2;
% Plot the tones' product and sum
figure(2)
subplot(2,1,1)
plot(n/Fs, Product, '-b'),
ylabel('Product'); xlabel('Time (sec.)');
grid on, zoom on
hold on
Red_Curve = 0.5*cos(2*pi*10*n/Fs) + 0.5; % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -1.25, 1.5]);
hold off
grid on, zoom on
subplot(2,1,2)
plot(n/Fs, Sum, '-b')
hold on
Red_Curve = 2*cos(2*pi*5*n/Fs); % Used for plotting only
plot(n/Fs, Red_Curve, '-r')
axis([0, 0.3, -2.4, 3]);
hold off
ylabel('Sum'); xlabel('Time (sec.)');
grid on, zoom on
% Play all the signals
sound(Wave_1, Fs)
pause(1.2)
sound(Wave_2, Fs)
pause(1.2)
sound(Product, Fs)
pause(1.2)
sound(Sum, Fs)
% Spec analysis of the "Sum" signal
Spec = fft(Sum);
Spec_Mag = abs(Spec);
Freq = n*Fs/N; % Freq vector in Hz
figure (3) % Plot positive-freq spec amg
plot(Freq(1:N/16), Spec_Mag(1:N/16))
title('Spec Mag of "Sum" signal')
ylabel('Mag.'), xlabel('Hz')
grid on, zoom on
```

## Computing the Nth Roots of a Number

```
function [nth_Roots] = roots_nth(Num, n, Plot)
% ROOTS_NTH computes the nth roots of Num.
%
% Inputs:
% Num is the number(s) for which the nth roots will be computed.
% Num can be real-valued or complex-valued.
% If Num has more than one element, then Num must be a row vector.
% n must be a positive integer.
% Plot is a string equal to 'Y' if a complex plane plot of
% the n complex roots of each element in Num is desired.
% Principle root is plotted as a blue square. Remaining roots
% are red dots. If input Plot string does not exist, then no
% plot is generated.
%
% Output:
% nth_Roots is a column vector (with n rows) giving the nth roots of Num.
%
% Example:
% clear, clc % clear all variables
% Num = [2+j*11, 3*exp(j*pi)]; % Two numbers
% n = 3; % find the 3rd roots
% [nth_Roots] = roots_nth(Num, n, 'Y')
%
% Results are:
% nth_Roots =
%
% 2.0000 + 1.0000i 0.7211 + 1.2490i
% -1.8660 + 1.2321i -1.4422 + 0.0000i
% -0.1340 - 2.2321i 0.7211 - 1.2490i
%
% [R. Lyons, Jan. 2013]
% Input argument check
if (nargin<2),
disp('Not enough input arguments!!')
return
end;
if (nargin==2)
Plot == 'N'
end
Num_of_Columns = max(size(Num)); % How many elements are in Num?
for Column_Loop = 1:Num_of_Columns;
Mag(Column_Loop) = abs(Num(Column_Loop));
Angle_Degrees(Column_Loop) = angle(Num(Column_Loop))*180/pi;
for k = 1:n
Root_Mag = Mag(Column_Loop)^(1/n);
Root_Angle_Deg(k) = Angle_Degrees(Column_Loop)/n + (k-1)*360/n;
nth_Roots(k,Column_Loop) = Root_Mag*exp(j*Root_Angle_Deg(k)*pi/180);
end
% Plot roots, on complex plane(s), if necessary
if Plot == 'Y'
figure
% Plot unit circle
Num_Pts = 200; % Number of circle points
Index = 0 : Num_Pts-1;
Alpha = 1/Num_Pts;
Points = Root_Mag*exp(j*2*pi*Alpha*Index);
plot(Points, ':k','markersize',5)
grid on
hold on
% Plot principle root (square blue dot) and list results in plot
axis([-1.1*Root_Mag, 1.1*Root_Mag, -1.1*Root_Mag, 1.1*Root_Mag]);
plot(Root_Mag*cos(Root_Angle_Deg(1)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(1)*pi/180), 'bs');
text(-Root_Mag/2, 3*Root_Mag/4, ['Magnitudes = ',num2str(Root_Mag)]);
text(-Root_Mag/2, 3*Root_Mag/4-Root_Mag/10, ...
['Principle Angle (deg.) = ',num2str(Root_Angle_Deg(1))]);
% Plot remaining roots as red dots
for k = 2:n
plot(Root_Mag*cos(Root_Angle_Deg(k)*pi/180), ...
Root_Mag*sin(Root_Angle_Deg(k)*pi/180), 'ro');
text(-Root_Mag/2, 3*Root_Mag/4-k*Root_Mag/10, ...
[' Angle ',num2str(k),' (deg.) = ', ...
num2str(Root_Angle_Deg(k))]);
end % End k loop
xlabel('Real'); ylabel('Imag.');
hold off
end % End 'if Plot == 'Y''
end % End Column_Loop
```

## Computing CIC Filter Register Pruning Using Matlab [Updated]

```
% Filename: CIC_Word_Truncation.m
%
% Computes CIC decimation filter accumulator register
% truncation in each filter stage based on Hogenauer's
% 'accumulator register pruning' technique.
%
% Inputs:
% N = number of decimation CIC filter stages (filter order).
% R = CIC filter rate change factor (decimation factor).
% M = differential delay.
% Bin = number of bits in an input data word.
% Bout = number of bits in the filter's final output data word.
% Outputs:
% Stage number (ranges from 1 -to- 2*N+1).
% Bj = number of least significant bits that can be truncated
% at the input of a filter stage.
% Accumulator widths = number of a stage's necessary accumulator
% bits accounting for truncation.
% Richard Lyons Feb., 2012
clear, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define CIC filter parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%N = 4; R = 25; M = 1; Bin = 16; Bout = 16; % Hogenauer paper, pp. 159
%N = 3; R = 32; M = 2; Bin = 8; Bout = 10; % Meyer Baese book, pp. 268
%N = 3; R = 16; M = 1; Bin = 16; Bout = 16; % Thorwartl's PDF file
N = 3; R = 8; M = 1; Bin = 12; Bout = 12; % Lyons' blog Figure 2 example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find h_sub_j and "F_sub_j" values for (N-1) cascaded integrators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
for j = N-1:-1:1
h_sub_j = [];
h_sub_j((R*M-1)*N + j -1 + 1) = 0;
for k = 0:(R*M-1)*N + j -1
for L = 0:floor(k/(R*M)) % Use uppercase "L" for loop variable
Change_to_Result = (-1)^L*nchoosek(N, L)...
*nchoosek(N-j+k-R*M*L,k-R*M*L);
h_sub_j(k+1) = h_sub_j(k+1) + Change_to_Result;
end % End "L" loop
end % End "k" loop
F_sub_j(j) = sqrt(sum(h_sub_j.^2));
end % End "j" loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for up to seven cascaded combs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j_for_many_combs = sqrt([2, 6, 20, 70, 252, 924, 3432]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for last integrator stage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(N) = F_sub_j_for_many_combs(N-1)*sqrt(R*M); % Last integrator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute F_sub_j for N cascaded filter's comb stages
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2*N:-1:N+1
F_sub_j(j) = F_sub_j_for_many_combs(2*N -j + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define "F_sub_j" values for the final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_sub_j(2*N+1) = 1; % Final output register truncation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of minus log base 2 of "F_sub_j" values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Minus_log2_of_F_sub_j = -log2(F_sub_j)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute total "Output_Truncation_Noise_Variance" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CIC_Filter_Gain = (R*M)^N;
Num_of_Bits_Growth = ceil(log2(CIC_Filter_Gain));
% The following is from Hogenauer's Eq. (11)
%Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin -1;
Num_Output_Bits_With_No_Truncation = Num_of_Bits_Growth +Bin;
Num_of_Output_Bits_Truncated = Num_Output_Bits_With_No_Truncation -Bout;
Output_Truncation_Noise_Variance = (2^Num_of_Output_Bits_Truncated)^2/12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute log base 2 of "Output_Truncation_Noise_Standard_Deviation" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output_Truncation_Noise_Standard_Deviation = ...
sqrt(Output_Truncation_Noise_Variance);
Log_base2_of_Output_Truncation_Noise_Standard_Deviation = ...
log2(Output_Truncation_Noise_Standard_Deviation);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute column vector of "half log base 2 of 6/N" terms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Half_Log_Base2_of_6_over_N = 0.5*log2(6/N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute desired "B_sub_j" vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B_sub_j = floor(Minus_log2_of_F_sub_j ...
+ Log_base2_of_Output_Truncation_Noise_Standard_Deviation ...
+ Half_Log_Base2_of_6_over_N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' '), disp(' ')
disp(['N = ',num2str(N),', R = ',num2str(R),', M = ',num2str(M),...
', Bin = ', num2str(Bin),', Bout = ',num2str(Bout)])
disp(['Num of Bits Growth Due To CIC Filter Gain = ', ...
num2str(Num_of_Bits_Growth)])
disp(['Num of Accumulator Bits With No Truncation = ', ...
num2str(Num_Output_Bits_With_No_Truncation)])
% disp(['Output Truncation Noise Variance = ', ...
% num2str(Output_Truncation_Noise_Variance)])
% disp(['Log Base2 of Output Truncation Noise Standard Deviation = ',...
% num2str(Log_base2_of_Output_Truncation_Noise_Standard_Deviation)])
% disp(['Half Log Base2 of 6/N = ', num2str(Half_Log_Base2_of_6_over_N)])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and display "Results" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for Stage = 1:2*N
Results(Stage,1) = Stage;
Results(Stage,2) = F_sub_j(Stage);
Results(Stage,3) = Minus_log2_of_F_sub_j(Stage);
Results(Stage,4) = B_sub_j(Stage);
Results(Stage,5) = Num_Output_Bits_With_No_Truncation -B_sub_j(Stage);
end
% Include final output stage truncation in "Results" matrix
Results(2*N+1,1) = 2*N+1; % Output stage number
Results(2*N+1,2) = 1;
Results(2*N+1,4) = Num_of_Output_Bits_Truncated;
Results(2*N+1,5) = Bout;
%Results % Display "Results" matrix in raw float-pt.form
% % Display "F_sub_j" values if you wish
% disp(' ')
% disp(' Stage Fj -log2(Fj) Bj Accum width')
% for Stage = 1:2*N+1
% disp([' ',sprintf('%2.2g',Results(Stage,1)),sprintf('\t'),sprintf('%12.3g',Results(Stage,2)),...
% sprintf('\t'),sprintf('%7.5g',Results(Stage,3)),sprintf('\t'),...
% sprintf('%7.5g',Results(Stage,4)),sprintf('\t'),sprintf('%7.5g',Results(Stage,5))])
% end
% Display Stage number, # of truncated input bits, & Accumulator word widths
disp(' ')
disp(' Stage(j) Bj Accum (adder) width')
for Stage = 1:2*N
disp([' ',sprintf('%2.0f',Results(Stage,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(Stage,4)),sprintf('\t'),...
sprintf('%7.5g',Results(Stage,5))])
end
disp([' ',sprintf('%2.0f',Results(2*N+1,1)),...
sprintf('\t'),...
sprintf('%5.5g',Results(2*N+1,4)),sprintf('\t'),...
sprintf('%7.5g',Results(2*N+1,5)),' (final truncation)'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```

## Adding a Controlled Amount of Noise to a Noise-Free Signal

```
function [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
% SNR_Set(x, Desired_SNR_dB) returns the real-valued
% input 'Signal' contaminated with normally-distributed,
% zero-mean, random noise. The signal-to-noise ratio
% (SNR in dB) of the output 'Noisy_Signal' signal is
% controlled by the input 'Desired_SNR_dB' argument measured
% in dB.
% Example:
% Npts = 128; % Number of time samples
% n = 0:Npts-1; % Time-domain index
% Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% Desired_SNR_dB = 3; % Set SNR of output 'Noisy_Signal' to +3 dB
% [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
% Author: Richard Lyons [December 2011]
%******************************************
Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one
Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));
K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10); % Scale factor
New_Noise = sqrt(K)*Noise; % Change Noise level
%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))
Noisy_Signal = Signal + New_Noise;
'SNR_Set()' Function Test Code:
% Filename SNR_Set_test.m
%
% Tests the 'SNR_Set()" function. Adds a predefined
% amount of random noise to a noise-free signal such that
% the noisy signal has a desired signal-to-noise ratio (SNR).
%
% Author: Richard Lyons [December 2011]
clear, clc
% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5; % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal
Desired_SNR_dB = 3 % Set desired SNR in dB
[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal
% Plot original and 'noisy' signals
figure(1)
subplot(2,1,1)
plot(n, Signal, '-bo', 'markersize', 2)
title('Original Signal')
grid on, zoom on
subplot(2,1,2)
plot(n, Noisy_Signal, '-bo', 'markersize', 2)
title('Noisy Signal')
xlabel('Time-samples')
grid on, zoom on
% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
figure(2)
plot(Spec_Mag, '-bo', 'markersize', 2)
title('Spec Mag of Noisy Signal')
xlabel('Freq-samples'), ylabel('Linear')
grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)
```

## Testing the Flat-Top Windowing Function

```
% Code for testing the 'Wind_Flattop(Spec)' function in
% reducing 'scalloping loss' errors in time signal amplitude
% estimation.
%
% Generates a time-domain sinusoid, computes its FFT, and passes
% that FFT sequence to the 'Wind_Flattop(Spec)' function.
%
% The maximum output sample of the 'Wind_Flattop(Spec)' function
% is used to estimate the peak amplitude of the original
% sinusoidal time-domain test signal.
%
% The user controls the values for the test sinusoid's
% 'Test_Freq' and peak amplitude 'Peak_Amp' near lines 23 & 24.
%
% Richard Lyons [December, 2011]
clear all, clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define test parameters
Test_Freq = 7.22; % Test tone's freq. Must be less that N/2
Peak_Amp = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 64; % Number of time samples
Index = (0:N-1);
X = Peak_Amp*cos(2*pi*(Test_Freq)*Index/N + pi/3);
figure(1), clf
subplot(2,1,1)
plot(X,':ko', 'markersize', 4)
title('Original time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT the input 'X' sequence and call 'Wind_Flattop()' function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Spec = fft(X);
[Windowed_Spec] = Wind_Flattop(Spec);
subplot(2,1,2), plot(abs(Spec),'ko', 'markersize', 3)
title('SpecMag of unwindowed time signal'), grid on, zoom on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results accuracy (error)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(['Test Freq = ',num2str(Test_Freq),...
', True Peak Amplitude = ',num2str(Peak_Amp)])
Mag_peak_unwindowed = max(abs(Spec));
Unwindowed_Amp_Estimate = 2*Mag_peak_unwindowed/N;
Unwindowed_Amp_Estimate_Error_in_dB = ...
20*log10(Unwindowed_Amp_Estimate/Peak_Amp);
disp(' ')
disp(['Unwindowed Peak Amplitude Estimate = ',...
num2str(Unwindowed_Amp_Estimate)])
disp(['Unwindowed Estimate Error in dB = ',...
num2str(Unwindowed_Amp_Estimate_Error_in_dB),' dB'])
M_peak_windowed = max(abs(Windowed_Spec));
Windowed_Amp_Estimated = 2*M_peak_windowed/N;
Windowed_Amp_Estimation_Error_in_dB = ...
20*log10(Windowed_Amp_Estimated/Peak_Amp);
disp(' ')
disp(['Windowed Peak Amplitude Estimate = ',...
num2str(Windowed_Amp_Estimated)])
disp(['Windowed Estimate Error in dB = ',...
num2str(Windowed_Amp_Estimation_Error_in_dB),' dB'])
```

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

```
function [Windowed_Spec] = Wind_Flattop(Spec)
% Given an input spectral sequence 'Spec', that is the
% FFT of some time sequence 'x', Wind_Flattop(Spec)
% returns a spectral sequence that is equivalent
% to the FFT of a flat-top windowed version of time
% sequence 'x'. The peak magnitude values of output
% sequence 'Windowed_Spec' can be used to accurately
% estimate the peak amplitudes of sinusoidal components
% in time sequence 'x'.
% Input: 'Spec' (a sequence of complex FFT sample values)
%
% Output: 'Windowed_Spec' (a sequence of complex flat-top
% windowed FFT sample values)
%
% Based on Lyons': "Reducing FFT Scalloping Loss Errors
% Without Multiplication", IEEE Signal Processing Magazine,
% DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
% Richard Lyons [December, 2011]
N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];
% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);
Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);
% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);
Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);
% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function
```

## Computing FFT Twiddle Factors

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

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

```
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
% by Jeff Miller. Downloaded from Matlab Central.
[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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
```