DSPRelated.com
   Code Snippets

Peak-to-average ratio analysis

Markus Nentwig ● 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

kaz - ● 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

Jesus Selva ● 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

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

Computing FFT Twiddle Factors

Rick Lyons ● November 28, 2011●10 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

Rick Lyons ● 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

Rick Lyons ● 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 
% 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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Jesus Selva ● 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

kaz - ● 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lsb = log2(M) - log2(k);                 %LSBs discarded when addressing
lut = round(A*sin(2*pi*(0:k-1)/k));      %lut, one cycle sine data

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

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

Complex filter with three multipliers per tap

kaz - ● 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
h1 = real(h)+imag(h);         %extra adder
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);