DSPRelated.com
   Code Snippets

Phase drift of NCO

kaz - ● December 10, 2011 ● Coded in Matlab 
clear all; close all;

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

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

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

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

Spectrum analyzer

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

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

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

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

    saPar = struct('rate_Hz', rate_Hz, ...
                   'pRef_W', 1e-3, ...
                   'fMin_Hz', 1e3, ...
                   'RBW_power_Hz', 1e5, ...
                   'RBW_window_Hz', 1e3, ...
                   'filter', 'none', ...
                   'plotStyle', 'b-', ...
                   'logscale', true, ...
                   'fig', 2, ...
                   'noisefloor_dB', -300);
    
    tmp1 = normalize(IIR_filterExample('signal', noise, fPar));
    tmp2 = normalize(IIR_filterExample('signal', pulse, fPar));
    
    spectrumAnalyzer('signal', tmp1, saPar);
    spectrumAnalyzer('signal', tmp2, saPar, 'plotStyle', 'k-');
end

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

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

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

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

    % clip to -1..1 range
    c=min(c, 1);
    c=max(c, -1);
    
    % raised-cosine filter
    H = 1/2+cos(pi/2*(c+1))/2;

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

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

% ************************************
% continuous-time filter example
% ************************************
function sig = IIR_filterExample(varargin)
    p = vararginToStruct(varargin);
    
    % design continuous-time IIR filter
    [IIR_b, IIR_a] = cheby1(p.chebyOrder, p.chebyRipple_dB, 1, 's');
    
    % evaluate on the FFT frequency grid
    fb_Hz = FFT_frequencyBasis(numel(p.signal), p.rate_Hz);
    
    % Laplace domain operator, normalized to filter cutoff frequency
    % (the above cheby1 call designs for unity cutoff)
    s = 1i * fb_Hz / p.fCorner_Hz;
    
    % polynomial in s
    H_num = polyval(IIR_b, s);
    H_denom = polyval(IIR_a, s);
    H = H_num ./ H_denom;
    
    % apply filter
    sig = ifft(fft(p.signal) .* H);
end

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

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

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

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

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

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

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

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

    winLen=0;
    if ~isempty(winLeft)  
        
        % zero-pad the power spectrum in the middle with a length
        % equivalent to the window size. 
        % this guarantees that there is enough bandwidth for the filter!
        % one window length is enough, the spillover from both sides overlaps
        % Scale accordingly.
        winLen=size(winLeft, 2)+size(winRight, 2);
        
        % note: fr is the power shown in the plot, NOT a frequency
        % domain representation of a signal.
        % There is no need to renormalize because of the length change
        center=floor(nSamples/2)+1;
        rem=nSamples-center;
        fr=[fr(1:center-1), zeros(1, winLen-1), fr(center:end)];
        % construct window with same length as fr
        win=zeros(size(fr));
        win(1:1+size(winLeft, 2)-1)=winLeft;
        win(end-size(winRight, 2)+1:end)=winRight;
        
        assert(isreal(win));
        assert(isreal(fr));
        assert(size(win, 2)==size(fr, 2));
        
        % convolve using FFT
        win=fft(win);
        fr=fft(fr);
        
        fr=fr .* win;
        fr=ifft(fr);
        fr=real(fr); % chop off roundoff error imaginary part
        fr=max(fr, 0); % chop off small negative numbers
        
        % remove padding
        fr=[fr(1:center-1), fr(end-rem+1:end)];
    end

    % ************************************
    % build frequency basis and rotate 0 Hz to center
    % ************************************
    fb = FFT_frequencyBasis(numel(fr), p.rate_Hz);    
    fr = fftshift(fr);
    fb = fftshift(fb);
    
    if false
        % use in special cases (very long signals)
        
        % ************************************
        % data reduction:
        % If a filter is used, details smaller than windowBW are
        % suppressed already by the filter, and using more samples
        % gives no additional information. 
        % ************************************
        if numel(fr) > p.NMax
            switch(p.filter)
              case 'gaussian'
                % 0.2 offset from the peak causes about 0.12 dB error
                incr=floor(windowBW/binspacing_Hz/5);
              case 'brickwall'
                % there is no error at all for a peak
                incr=floor(windowBW/binspacing_Hz/3);
              case 'none'
                % there is no filter, we cannot discard data at this stage
                incr=-1;
            end
            
            if incr > 1 
                fr=fr(1:incr:end);
                fb=fb(1:incr:end);
            end  
        end
    end
    
    % ************************************
    % data reduction: 
    % discard beyond fMin / fMax, if given
    % ************************************
    indexMin = 1;
    indexMax = numel(fb);  
    flag=0;
    if isfield(p, 'fMin_Hz')
        indexMin=min(find(fb >= p.fMin_Hz));
        flag=1;
    end
    if isfield(p, 'fMax_Hz')
        indexMax=max(find(fb <= p.fMax_Hz));
        flag=1;
    end
    if flag
        fb=fb(indexMin:indexMax);
        fr=fr(indexMin:indexMax);
    end
    
    if p.NMax > 0
        if p.logscale
            % ************************************
            % Sample number reduction for logarithmic
            % frequency scale
            % ************************************     
            assert(isfield(p, 'fMin_Hz'), 'need fMin_Hz in logscale mode');
            assert(p.fMin_Hz > 0, 'need fMin_Hz > 0 in logscale mode');
            if ~isfield(p, 'fMax_Hz')
                p.fMax_Hz = p.rate_Hz / 2;
            end
            
            % averaged output arrays
            fbOut=zeros(1, p.NMax-1);
            frOut=zeros(1, p.NMax-1);

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

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

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

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

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

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

    % ************************************
    % Post-processing: 
    % Translate frequency axis to RF
    % ************************************
    fb = fb + p.originMapsTo_Hz;
    
    % ************************************
    % convert to requested units
    % ************************************
    switch(p.freqUnit)
      case 'Hz'
      case 'kHz'
        fb = fb / 1e3;
      case 'MHz'
        fb = fb / 1e6;
      case 'GHz'
        fb = fb / 1e9;
      otherwise
        error('unsupported frequency unit');
    end
    
    % ************************************
    % Plot (if requested)
    % ************************************    
    if p.fig > 0
        % *************************************************************
        % title
        % *************************************************************
        if isfield(p, 'title')
            t=['"', p.title, '";'];
        else
            t='';
        end
        switch(p.filter)
          case 'gaussian'
            t=[t, ' filter: Gaussian '];
          case 'brickwall'
            t=[t, ' filter: ideal bandpass '];
          case 'none'
            t=[t, ' filter: none '];
          otherwise
            assert(0)
        end
        
        if ~strcmp(p.filter, 'none')
            t=[t, '(', mat2str(p.RBW_window_Hz), ' Hz)'];       
        end
        
        if isfield(p, 'pNom_dBm')
            % *************************************************************
            % show dB relative to a given absolute power in dBm 
            % *************************************************************
            fr=fr-p.pNom_dBm;
            yUnit='dB';
        else
            yUnit='dBm';
        end
        
        % *************************************************************
        % plot
        % *************************************************************
        figure(p.fig);    
        if p.logscale
            handle = semilogx(fb, fr, p.plotStyle);
        else
            handle = plot(fb, fr, p.plotStyle);
        end
        hold on; % after plot, otherwise prevents logscale
        
        if isfield(p, 'lineWidth')
            set(handle, 'lineWidth', p.lineWidth);
        end
        
        % *************************************************************
        % axis labels
        % *************************************************************
        xlabel(p.freqUnit);
        if singleBinMode
            ylabel([yUnit, ' (continuous wave)']);    
        else
            ylabel([yUnit, ' in ', mat2str(p.RBW_power_Hz), ' Hz integration BW']);
        end
        title(t);
        
        % *************************************************************
        % adapt y range, if requested
        % *************************************************************
        y=ylim();
        y1=y(1); y2=y(2);

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

function [winLeft, winRight] = spectrumAnalyzerGetWindow(filter, singleBinMode, RBW_window_Hz, binspacing_Hz)
    
    switch(filter)
      case 'gaussian'
        
        % construct Gaussian filter
        % -60 / -3 dB shape factor 4.472
        nRBW=6;
        nOneSide=ceil(RBW_window_Hz/binspacing_Hz*nRBW); 
        
        filterBase=linspace(0, nRBW, nOneSide);
        winLeft=exp(-(filterBase*0.831) .^2);
        winRight=fliplr(winLeft(2:end)); % omit 0 Hz bin
        
      case 'brickwall'
        nRBW=1;
        n=ceil(RBW_window_Hz/binspacing_Hz*nRBW); 
        n1 = floor(n/2);
        n2 = n - n1;
        
        winLeft=ones(1, n1);
        winRight=ones(1, n2); 
        
      case 'none'
        winLeft=[];
        winRight=[];
      
      otherwise
        error('unknown RBW filter type');
    end
    
    % the window is not supposed to shift the spectrum.
    % Therefore, the first bin is always in winLeft (0 Hz):
    if size(winLeft, 2)==1 && isempty(winRight)
        % there is no use to convolve with one-sample window
        % it's always unity
        winLeft=[];
        winRight=[];
        tmpwin=[];
    end
    
    if ~isempty(winLeft)
        % (note: it is not possible that winRight is empty, while winLeft is not)
        if singleBinMode
            % normalize to unity at 0 Hz
            s=winLeft(1);
        else
            % normalize to unity area under the filter
            s=sum(winLeft)+sum(winRight);
        end
        winLeft=winLeft / s;
        winRight=winRight / s;
    end
end

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

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

Peak-to-average ratio analysis

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

Tremolo audio effect (amplitude modulation)

Gabriel Rivas ● November 29, 2011●3 comments ● Coded in C 
/************************Tremolo.c*************************/

#include "Tremolo.h"

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

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

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

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

void Tremolo_sweep(void) {

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

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

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

#endif

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

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

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

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

Auto Wah and Envelope follower audio effect

Gabriel Rivas ● November 29, 2011 ● Coded in C 
#include "bp_iir.h"
#include "AutoWah.h"

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

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

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

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

    bp_iir_init(sampling,gainfactor,Q,freq_step,minf);
	f_step = freq_step; 
        
	//Lowpass filter parameters	
	/*
	   C = 1/TAN(PI*Fc/Fs)
	*/
	C = 1018.59;

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

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

double AutoWah_process(double xin) {
	double yout;

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

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

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

	if (!--counter) {
		/*The output of the LPF (y[2]) is scaled by 0.1
		in order to be used as a LFO to control the band pass filter
		*/
	    filter_index = (double)min_freq + y[2]*0.1;
                        
            /*The scaling value is determined as a value
            that would keep the filter index within the
            range of max_freq
	    */
            if(filter_index > max_freq) {
                filter_index = max_freq;
            } 
            if(filter_index < 0) {
                filter_index = 0;
            }

	    bp_iir_setup(&H,filter_index);

	    counter = counter_limit; 
	}
}

Computing FFT Twiddle Factors

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