Code Snippets

Beat Notes

Rick Lyons ● March 13, 2013 ● Coded in Matlab 
%     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

Rick Lyons ● January 9, 2013●1 comment ● Coded in Matlab 
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

Delay estimation revisited

Markus Nentwig ● June 9, 2012●5 comments ● Coded in Matlab 
% ****************************************************************
% find sub-sample delay and scaling factor between two cyclic signals 
% to maximize crosscorrelation
% Markus Nentwig, 120609_v1.1
% ****************************************************************
function iterDelayEstDemo();
    close all;
    n = 1024;    
    % n = 1024 * 256; disp('*** test: long signal enabled ***');

    % ****************************************************************
    % random signal
    % ****************************************************************
    fd = randn(1, n) + 1i * randn(1, n);
    
    % ****************************************************************
    % lowpass filter
    % ****************************************************************
    f = binFreq(n);
    fd(abs(f) > 0.045) = 0;
    s1 = real(ifft(fd)) * sqrt(n);
    
    % ****************************************************************
    % create delayed 2nd signal
    % ****************************************************************
    dTest_samples = 12.3456;
    cTest = 1.23456;
    % cTest = cTest + 1i; disp('*** test: complex coeff enabled ***'); 
    % cTest = -cTest; disp('*** test: negative coeff enabled ***'); 

    s2 = cTest * cs_delay(s1, 1, dTest_samples);
    % s2 = s2 + 0.5*randn(size(s2)); disp('*** test: noise enabled ***');
    
    % ****************************************************************
    % estimate delay
    % ****************************************************************
    [delay_samples, coeff] = iterDelayEst(s1, s2);

    % ****************************************************************
    % correct it
    % ****************************************************************
    s2a = cs_delay(s1, 1, delay_samples);
    s2b = s2a * coeff;

    figure(); hold on;
    h = plot(real(s1), 'k'); set(h, 'lineWidth', 3);
    h = plot(real(s2), 'b'); set(h, 'lineWidth', 3);
    h = plot(real(s2a), 'r'); set(h, 'lineWidth', 1);
    h = plot(real(s2b), 'm'); set(h, 'lineWidth', 1);
    xlim([1, numel(s1)]);
    xlabel('samples');
    legend('s1', 's2', 's2 un-delayed', 's2 un-delayed and scaled');
    title('test signals');
    
    format long;
    disp('nominal delay of s2 relative to s1')';
    dTest_samples
    disp('iterDelayEst() returned:');
    delay_samples
    disp('original scaling factor:');
    cTest
    disp('estimated scaling factor:');
    coeff
end

% ****************************************************************
% estimates delay and scaling factor 
% ****************************************************************
function [delay_samples, coeff] = iterDelayEst(s1, s2)
    
    s1 = s1(:) .'; % force row vectors
    s2 = s2(:) .';
    rflag = isreal(s1) && isreal(s2);
    n = numel(s1);
    halfN = floor(n/2);
    assert(numel(s2) == n, 'signals must have same length');

    % ****************************************************************
    % constants
    % ****************************************************************    
    % exit if uncertainty below threshold
    thr_samples = 1e-7;

    % exit after fixed number of iterations
    nIter = 25;

    % frequency domain representation of signals
    fd1 = fft(s1);
    fd2 = fft(s2);    

    % first round: No delay was applied
    tau = [];
    fd2Tau = fd2; % delayed s2 in freq. domain
    
    % frequency corresponding to each FFT bin -0.5..0.5
    f = binFreq(n);

    % uncertainty plot data
    e = [];

    % normalization factor
    nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n; % normalizes to 1
    
    % search window: 
    % known maximum and two surrounding points
    x1 = -1;
    x2 = -1;
    x3 = -1;
    y1 = -1;
    y2 = -1;
    y3 = -1;
    
    % ****************************************************************
    % iteration loop
    % ****************************************************************
    for count = 1:nIter
    
        % ****************************************************************
        % crosscorrelation with time-shifted signal
        % ****************************************************************
        xcorr = abs(ifft(fd2Tau .* conj(fd1)))/ nf;

        % ****************************************************************
        % detect peak
        % ****************************************************************
        if isempty(tau)

            % ****************************************************************
            % startup
            % initialize with three adjacent bins around peak
            % ****************************************************************
            ix = find(xcorr == max(xcorr));
            ix = ix(1); % use any, if multiple bitwise identical peaks
            
            % indices of three bins around peak
            ixLow = mod(ix-1-1, n) + 1; % one below
            ixMid = ix;
            ixHigh = mod(ix-1+1, n) + 1; % one above

            % delay corresponding to the three bins
            tauLow = mod(ixLow -1 + halfN, n) - halfN;
            tauMid = mod(ixMid -1 + halfN, n) - halfN;         
            tauHigh = mod(ixHigh -1 + halfN, n) - halfN; 

            % crosscorrelation value for the three bins
            xcLow = xcorr(ixLow);
            xcMid = xcorr(ixMid);
            xcHigh = xcorr(ixHigh);
            
            x1 = tauLow;
            x2 = tauMid;
            x3 = tauHigh;
            y1 = xcLow;
            y2 = xcMid; 
            y3 = xcHigh;
        else
            % ****************************************************************
            % only main peak at first bin is of interest
            % ****************************************************************
            tauMid = tau;
            xcMid = xcorr(1);

            if xcMid > y2
                % ****************************************************************
                % improve midpoint
                % ****************************************************************
                if tauMid > x2
                    % midpoint becomes lower point
                    x1 = x2;
                    y1 = y2;
                else
                    % midpoint becomes upper point
                    x3 = x2;
                    y3 = y2;
                end
                x2 = tauMid;
                y2 = xcMid;
            
            elseif tauMid < x2
                % ****************************************************************
                % improve low point
                % ****************************************************************
                assert(tauMid >= x1); % bitwise identical is OK
                assert(tauMid > x1 || xcMid > y1); % expect improvement
                x1 = tauMid;
                y1 = xcMid;
            elseif tauMid > x2 
                % ****************************************************************
                % improve high point
                % ****************************************************************
                assert(tauMid <= x3); % bitwise identical is OK                
                assert((tauMid < x3) || (xcMid > y3)); % expect improvement
                x3 = tauMid;
                y3 = xcMid;
            else
                assert(false, '?? evaluated for existing tau ??');
            end
        end

        % ****************************************************************
        % calculate uncertainty (window width)
        % ****************************************************************
        eIter = abs(x3 - x1);
        e = [e, eIter];
        if eIter < thr_samples
            % disp('threshold reached, exiting');
            break;
        end

        if y1 == y2 || y2 == y3
            % reached limit of numerical accuracy on one side
            usePoly = 0;
        else
            % ****************************************************************
            % fit 2nd order polynomial and find maximum
            % ****************************************************************
            num = (x2^2-x1^2)*y3+(x1^2-x3^2)*y2+(x3^2-x2^2)*y1;
            denom = (2*x2-2*x1)*y3+(2*x1-2*x3)*y2+(2*x3-2*x2)*y1;
            if denom ~= 0
                tau = num / denom;
                % is the point within [x1, x3]? 
                usePoly = ((tau > x1) && (tau < x3));
            else
                usePoly = 0;
            end            
        end
        if ~usePoly
            % revert to linear interpolation on the side with the
            % less-accurate outer sample 
            % Note: There is no guarantee that the side with the more accurate
            % outer sample is the right one, as the samples aren't 
            % placed on a regular grid!
            % Therefore, iterate to improve the "worse" side, which will
            % eventually become the "better side", and iteration converges.
            
            tauLow = (x1 + x2) / 2;
            tauHigh = (x2 + x3) / 2;
            if y1 < y3
                o = [tauLow, tauHigh];
            else
                o = [tauHigh, tauLow];                
            end
            % don't choose point that is identical to one that is already known
            tau = o(1);
            if tau == x1 || tau == x2 || tau == x3
                tau = o(2);
                if tau == x1 || tau == x2 || tau == x3
                    break;
                end
            end
        end

        % ****************************************************************
        % advance 2nd signal according to location of maximum
        % phase shift in frequency domain - delay in time domain
        % ****************************************************************
        fd2Tau = fd2 .* exp(2i * pi * f * tau);
    end % for

    % ****************************************************************
    % plot the uncertainty (window size) over the number of iterations
    % ****************************************************************
    if true
        figure(); semilogy(e, '+-'); grid on;
        xlabel('iteration');
        title('uncertainty in delay');
    end

    % ****************************************************************
    % the delay estimate is the final location of the delay that 
    % maximized crosscorrelation (center of window).
    % ****************************************************************
    delay_samples = x2;

    % ****************************************************************
    % Coefficient: Turn signal 1 into signal 2
    % ****************************************************************
    coeff = fd2Tau * fd1' ./ (fd1 * fd1');

    % ****************************************************************
    % chop roundoff error, if input signals are known to be 
    % real-valued.
    % ****************************************************************
    if rflag
        coeff = real(coeff);
    end
end

% ****************************************************************
% frequency corresponding to FFT bin
% ****************************************************************
function f = binFreq(n)
    f = (mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
end

% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
    rflag = isreal(waveform);
    
    n = numel(waveform);
    cycLen_s = n / rate_Hz;
    nCyc = delay_s / cycLen_s();
    
    f = 0:(n - 1);
    f = f + floor(n / 2);
    f = mod(f, n);
    f = f - floor(n / 2);
    
    phase = -2 * pi * f * nCyc;
    rot = exp(1i*phase);
    
    waveform = ifft(fft(waveform) .* rot);
    if rflag
        waveform = real(waveform);
    end
end

Circular Convolution

Kadhiem Ayob ● March 3, 2012 ● Coded in Matlab 
%circular convolution
%for testing you may use:
h = fir1(20,.3);
x = randn(1,1024);

%function y = conv_circ(h,x)

y = conv(h,x);

L1 = length(h);
L2 = length(x);

%add end to start, add start to end
temp = y(1:L1-1);
y(1:L1-1) = y(1:L1-1) + y(L2+(1:L1-1));
y(L2+(1:L1-1)) = y(L2+(1:L1-1)) + temp;

%compare to direct convolution
y2 = conv(h,x);
plot(y,'o-');hold
plot(y2,'r.--')
legend('circular','direct')

Power Computation of a Digital Stream

Kadhiem Ayob ● February 26, 2012 ● Coded in Matlab 
clear all;

%example vector having ac & dc power
x = complex(randn(1,2048)+.2,randn(1,2048)-.1);

%total power, three equivalent methods
pwr1_total = mean(abs(x).^2);                %mean of squared values
pwr2_total = mean(real(x).^2 + imag(x).^2);
pwr3_total = mean(x .* conj(x));

%total expressed as rms
rms_total = sqrt(pwr1_total);

%dc power
pwr1_dc = mean(real(x))^2 + mean(imag(x))^2; %square of mean of values

%ac power
pwr1_ac = pwr1_total - pwr1_dc;              %mean of squares - square of mean

%relation of ac power to statistical variance 
pwr2_ac = var(x);              %approximately

%ac expressed as rms
rms1_ac = sqrt(pwr1_ac);

%ac relation to standard of deviation, std = sqrt(var)
rms2_ac = std(x);              %approximately

%dc relation to variance
pwr2_dc = pwr1_total - var(x); %approximately

fprintf('----------------------------------------------------\r');
fprintf('power(total),          : %1.5f, %1.5f, %1.5f\r',pwr1_total,pwr2_total,pwr3_total);
fprintf('rms(total)             : %1.5f\r',rms_total);
fprintf('power(ac),  variance   : %1.5f, %1.5f\r',pwr1_ac,pwr2_ac);
fprintf('rms(ac),    std        : %1.5f, %1.5f\r',rms1_ac,rms2_ac);
fprintf('power(dc), (total-var) : %1.5f, %1.5f\r',pwr1_dc,pwr2_dc);

Computing CIC Filter Register Pruning Using Matlab

Rick Lyons ● February 15, 2012●8 comments ● Coded in Matlab 
%  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.0g',Results(Stage,1)),...
        sprintf('\t'),...
        sprintf('%5.5g',Results(Stage,4)),sprintf('\t'),...
        sprintf('%7.5g',Results(Stage,5))])
end
    disp(['  ',sprintf('%2.0g',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)'])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Reading text files in Matlab

Kadhiem Ayob ● February 13, 2012 ● Coded in Matlab 
clear all;

% create test file
data = [1 5 2 4 3 3 4 2 5 1];
filename = 'test_file.txt';
fid = fopen(filename, 'w');
fprintf(fid, '%d %d\n', data);
fclose(fid);

% (1) read file using fscanf 
fid = fopen(filename, 'r'); 
y1 = fscanf(fid, '%d\n'); %interleaves columns
fclose(fid);

% (2) read file using textread (or textscan)
[ya,yb] = textread(filename,'%d%d');
y2 = [ya yb];

% (3) read file using importdata 
y3 = importdata(filename);

% (4) read file using load 
y4 = load(filename);

disp('-------------------------------------------------------------')
disp('    original vector data')
disp(data)
disp('    file content using fprintf')
disp(y2)
disp('    vector created by fscanf')
disp(y1)
disp('     matrix created by:')
disp('     textread   importdata    load')
disp([y2 y3 y4])

Ideal interpolation filter

Kadhiem Ayob ● January 21, 2012 ● Coded in Matlab 
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,1024));

up = 3;             %Interpolation factor
cutoff = .3;

intorder = 6;
h1 = intfilt(up, intorder, cutoff); %ideal filter
h1 = up*h1/sum(h1);

h2 = fir1(2*up*intorder-2,cutoff);  %ordinary LPF 
h2 = up*h2/sum(h2);

%upsample
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;

x_f1 = filter(h1,1,x_up);
x_f2 = filter(h2,1,x_up);

figure;
subplot(3,1,1);hold
plot(x_f1,'o--');
plot(x_f2,'r.-');
legend('ideal output','fir1 output');

subplot(3,1,2);
plot(x_f1-x_f2);
legend('error');

subplot(3,1,3);hold
plot(h1,'.-');
plot(h2,'r.-');
legend('ideal filter','fir1 filter');

checking resampling in time domain

Kadhiem Ayob ● January 21, 2012 ● Coded in Matlab 
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,512));
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor

%%%%%%%%%%%%%%%%%%%%% up/dn model %%%%%%%%%%%%%%%%%%%%%%
%upsample input
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f = filter(up*h,1,x_up);

%downsample signal by decimation
x_dn = x_f(1:dn:end);

delay = 30/2+1;  
figure;
subplot(2,1,1);hold;
plot(x_up,'o--');
plot(x_f(delay:end),'r.-');
legend('input with zeros','upsampled stage');

subplot(2,1,2);hold;
plot(x_up(1:dn:end),'o--');
plot(x_dn(ceil(delay/dn):end),'r.-');
legend('input with zeros','final signal');

Resampling filter performance

Kadhiem Ayob ● January 14, 2012●3 comments ● Coded in Matlab 
%%%%%%%%%%%%%%%%%% inputs for model %%%%%%%%%%%%%%%%
clear all; close all;

%example bandlimited impulse input & parameters
x = filter(fir1(70,.1),1,[1 zeros(1,2^15-1)]);
Fs = 120;           %MHz original sample rate             
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor
Fc = 12;             %MHz band centre (-Fs/2 ~ +Fs/2)
Fch = 0;             %MHz filter centre (-Fs*up/2 ~ +Fs*up/2)

%move signal to its centre
x = x.*exp(j*2*pi*(0:length(x)-1)*Fc/Fs);

%shift filter 
h = h.*exp(j*2*pi*(0:length(h)-1)*Fch/(Fs*up));

%%%%%%%%%%%%%%%%%%%%% model %%%%%%%%%%%%%%%%%%%%%%
%check signal in upsampled domain
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
[P, F] = pwelch(x_up, [], 0, 2^16, Fs*up,'twosided'); 
F = F - max(F)/2; 
P = fftshift(P);
y(find(P == 0)) = -100;  %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P));  %normalise

%check filter response in upsampled domain
H = fftshift(20*log10(abs(fft(h,2^16))));

subplot(2,1,1);
hold;grid;
plot(F, P_dB,'.-'); 
plot(F,H,'m--');
axis([min(F)-1 max(F)+1 -80 1]);
legend('upsampled signal','upsampling filter');

%check signal in downsampled domain
x_f = filter(h,1,x_up);
x_dn = x_f(1:dn:end);
[P, F] = pwelch(x_dn, [], 0, 2^16, Fs*up/dn,'twosided'); 
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100;   %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise

subplot(2,1,2)
plot(F,P_dB,'r.-');
grid;
axis([min(F)-1 max(F)+1 -80 1]);
legend('downsampled signal');