DSPRelated.com
   Code Snippets

Adaptive IIR filter using fixed point integer arithmetic

● March 23, 2013 ● Coded in C 
typedef signed long s32;

#define ABS(x)                 ((x) < 0 ? -(x) : (x))
#define CLAMP(x,min,max)       {if ((x) <= (min)) (x) = (min); else if ((x) >= (max)) (x) = (max);}

#define KLPF_MAX_ADAPT			232L
#define KLPF_MIN_ADAPT			 0L
#define KLPF_DEFAULT_ADAPT		160L

#define ACCEL_MIN_ADAPT     5
#define ACCEL_MAX_ADAPT     15

s32 gnKLpf          = KLPF_DEFAULT_ADAPT;   // IIR filter coefficient, 0 to 255
s32 gnKLpfMax       = KLPF_MAX_ADAPT;  
s32 gnAccel = 0;   // rate of change of climbrate
s32 gnCpsx256 = 0; // climbrate in centimeters per second, with 8 bit fractional precision
s32 gnCps; // climbrate in centimeters per second

// IIR low pass filter using fixed point integer arithmetic. 8bits of fractional precision for IIR coefficient K.
// So the actual floating point coefficient is k = K/256, 0.0 <= k < 1.0. The filter computes
// climbRateFiltered = climbRateFiltered*k + newClimbRate*(1.0 - k). So K values close to 256 will result in 
// heavy damping, K values close to 0 will result in almost no damping.  The IIR low pass filtered output gnCps 
// is the rounded up integer value, the 8bit fraction precision value gnCpsx256 is a global variable, so is
// maintained between calls to the filter. No integer divides required. 

void sns_LPFClimbRate(void)   {
   	s32 tmp;
   	tmp = gnCpsx256*gnKLpf + (gnCps<<8)*(256L-gnKLpf);
   	gnCpsx256 = (tmp >= 0 ? ((tmp + 128L)>>8) : ((tmp - 128L)>>8));
   	gnCps = (gnCpsx256>>8);
   	}  

// Adapt the IIR filter coeffient to the rate of change. If samples are not changing much, increase the damping.
// If the samples are changing by a large amount, reduce the damping. This is done to reduce the response delay to a
// a step change, while still being able to heavily damp out jitter in steady state noisy samples.
// This function basically linearly interpolates KLpf between KLPF_MAX and KLPF_MIN based on the acceleration.
// Values of acceleration outside experimentally determined limits are clamped to the limits, depends on the 
// application. A variable is used for KLpfMax to allow user to adjust the maximum allowed damping.

s32 sns_AdaptKLpf(void) {
	s32 klpf,nAccel;
	nAccel = ABS(gnAccel);
	CLAMP(nAccel, ACCEL_MIN_ADAPT, ACCEL_MAX_ADAPT);
	klpf = gnKLpfMax + ((KLPF_MIN_ADAPT-gnKLpfMax)*(nAccel - ACCEL_MIN_ADAPT))/(ACCEL_MAX_ADAPT-ACCEL_MIN_ADAPT);
	CLAMP(klpf, KLPF_MIN_ADAPT, KLPF_MAX_ADAPT);
	return klpf;
	}

// usage : 
// For each new data sample that arrives
//   1. calculate new unfiltered climbrate nCps from sample data buffer	
//	 2. gnAccel = (nCps - gnCps);
//	 3. (optionally low pass filter gnAccel using same type of IIR filter, if required)
//   4. gnCps = nCps;
//	 5. gnKLpf = sns_AdaptKLpf();
//	 6. sns_LPFClimbRate();

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

Power spectra of MPSK

Senthilkumar ● March 28, 2012●1 comment ● Coded in Scilab 
function[SB_PSK] = PowerSpectra_MPSK()
rb = input('Enter the bit rate=');
Eb = input('Enter the energy of the bit=');
f = 0:1/100:rb;
Tb = 1/rb;  //Bit duration
M = [2,4,8];
for j = 1:length(M)
  for i= 1:length(f)
    SB_PSK(j,i)=2*Eb*(sinc_new(f(i)*Tb*log2(M(j)))^2)*log2(M(j));
  end
end
a=gca();
plot2d(f*Tb,SB_PSK(1,:)/(2*Eb))
plot2d(f*Tb,SB_PSK(2,:)/(2*Eb),2)
plot2d(f*Tb,SB_PSK(3,:)/(2*Eb),5)
xlabel('Normalized Frequency ---->')
ylabel('Normalized Power Spectral Density--->')
title('Power Spectra of M-ary signals for M =2,4,8')
legend(['M=2','M=4','M=8'])
xgrid(1) 
endfunction
//Result
//Enter the bit rate in bits per second:2
//Enter the Energy of bit:1

Power spectra of MSK & QPSk

Senthilkumar ● March 28, 2012●2 comments ● Coded in Scilab 
function [SB_MSK,SB_QPSK]= PowerSpectra_MSK_QPSK()
//Comparison of QPSK and MSK Power Spectrums
rb = input('Enter the bit rate in bits per second:');
Eb = input('Enter the Energy of bit:');
f = 0:1/(100*rb):(4/rb);
Tb = 1/rb; //bit duration in seconds
for i = 1:length(f)
  if(f(i)==0.5)
    SB_MSK(i) = 4*Eb*f(i);
  else
    SB_MSK(i) = (32*Eb/(%pi^2))*(cos(2*%pi*Tb*f(i))/((4*Tb*f(i))^2-1))^2;
  end
    SB_QPSK(i)= 4*Eb*sinc_new((2*Tb*f(i)))^2;
end
a = gca();
plot(f*Tb,SB_MSK/(4*Eb));
plot(f*Tb,SB_QPSK/(4*Eb));
poly1= a.children(1).children(1);
poly1.foreground = 3;
xlabel('Normalized Frequency ---->')
ylabel('Normalized Power Spectral Density--->')
title('QPSK Vs MSK Power Spectra Comparison')
legend(['Minimum Shift Keying','QPSK'])
xgrid(1)
endfunction
//Result
//Enter the bit rate in bits per second:2
//Enter the Energy of bit:1

Power Spectrum of Discrete PAM signals

Senthilkumar ● March 28, 2012 ● Coded in Scilab 
function [Sxxf_NRZ_P,Sxxf_NRZ_BP,Sxxf_NRZ_UP,Sxxf_Manch]=PowerSpectra_PAM()
a = input('Enter the Amplitude value:');
fb = input('Enter the bit rate:');
Tb = 1/fb;  //bit duration
f = 0:1/(100*Tb):2/Tb;
for i = 1:length(f)
  Sxxf_NRZ_P(i) = (a^2)*Tb*(sinc_new(f(i)*Tb)^2);
  Sxxf_NRZ_BP(i) = (a^2)*Tb*((sinc_new(f(i)*Tb))^2)*((sin(%pi*f(i)*Tb))^2);
  if (i==1)
    Sxxf_NRZ_UP(i) = (a^2)*(Tb/4)*((sinc_new(f(i)*Tb))^2)+(a^2)/4;
  else
    Sxxf_NRZ_UP(i) = (a^2)*(Tb/4)*((sinc_new(f(i)*Tb))^2);
  end
  Sxxf_Manch(i) = (a^2)*Tb*(sinc_new(f(i)*Tb/2)^2)*(sin(%pi*f(i)*Tb/2)^2);
end
//Plotting
a = gca();
plot2d(f,Sxxf_NRZ_P)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
plot2d(f,Sxxf_NRZ_BP,2)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
plot2d(f,Sxxf_NRZ_UP,5)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
plot2d(f,Sxxf_Manch,9)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
xlabel('f*Tb------->')
ylabel('Sxx(f)------->')
title('Power Spectral Densities of Different Line Codinig Techniques')
xgrid(1)
legend(['NRZ Polar Format','NRZ Bipolar format','NRZ Unipolar format','Manchester format']);
endfunction
//Result
//Enter the Amplitude value:1
//Enter the bit rate:1

Raised Cosine Spectrum

Senthilkumar ● March 28, 2012 ● Coded in Scilab 
function [p]= RaisedCosineSpectrum()
//Practical Solution for Intersymbol Interference
//Raised Cosine Spectrum
rb = input('Enter the bit rate:');
Tb =1/rb;
t =-3:1/100:3;
Bo = rb/2;
Alpha =0;      //roll-off factor Intialized to zero
x =t/Tb;
for j =1:3
  for i =1:length(t)
    if((j==3)&((t(i)==0.5)|(t(i)==-0.5)))
        p(j,i) = sinc_new(2*Bo*t(i));
    else
        num =  sinc_new(2*Bo*t(i))*cos(2*%pi*Alpha*Bo*t(i));
        den =   1-16*(Alpha^2)*(Bo^2)*(t(i)^2)+0.01; 
        p(j,i)= num/den;
    end
  end
  Alpha = Alpha+0.5;
end
a =gca();
plot2d(t,p(1,:))
plot2d(t,p(2,:))
poly1= a.children(1).children(1);
poly1.foreground=2;
plot2d(t,p(3,:))
poly2= a.children(1).children(1);
po1y2.foreground=4;
poly2.line_style = 3;
xlabel('t/Tb------>');
ylabel('p(t)------->');
title('RAISED COSINE SPECTRUM - Practical Solution for ISI')
legend(['ROlloff Factor =0','ROlloff Factor =0.5','ROlloff Factor =1'])
xgrid(1)
endfunction
//Result
//Enter the bit rate:1

Modified Duobinary Signaling-Spectrum

Senthilkumar ● March 28, 2012 ● Coded in Scilab 
function[Amplitude_Response,Phase_Response]=Modified_Duobinary_Signaling()
//Modified Duobinary Signaling Scheme
//Magnitude and Phase Response
rb =  input('Enter the bit rate=');
Tb =1/rb;  //Bit duration
f = -rb/2:1/100:rb/2;
Amplitude_Response = abs(2*sin(2*%pi*f.*Tb));
Phase_Response = -(2*%pi*f.*Tb);
subplot(2,1,1)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot2d(f,Amplitude_Response,2)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
xlabel('Frequency f---->')
ylabel('|H(f)| ----->')
title('Amplitude Repsonse of Modified Duobinary Singaling')
xgrid(1)
subplot(2,1,2)
a=gca();
a.x_location ="origin";
a.y_location ="origin";
plot2d(f,Phase_Response,5)
poly1= a.children(1).children(1);  
poly1.thickness = 2;  // the tickness of a curve.
xlabel('                                           Frequency f---->')
ylabel('                                            <H(f) ----->')
title('Phase Repsonse of Modified Duobinary Singaling')
xgrid(1)
//Result
//Enter the bit rate=8 
endfunction
//Result
//Enter the bit rate=8

Duobinary Encoder and Decoder

Senthilkumar ● March 28, 2012 ● Coded in Scilab 
function [c,b_r]=Duobinary_EncDec(b)

// Precoded Duobinary coder and decoder
//b = input binary sequence:precoder input
//c = duobinary coder output
//b_r = duobinary decoder output
a(1) = xor(1,b(1));
if(a(1)==1)
  a_volts(1) = 1;
end
for k =2:length(b)
  a(k) = xor(a(k-1),b(k));
  if(a(k)==1)
    a_volts(k)=1;
  else
    a_volts(k)=-1;
  end
end
a = a';
a_volts = a_volts';
disp(a,'Precoder output in binary form:')
disp(a_volts,'Precoder output in volts:')
//Duobinary coder output in volts
c(1) = 1+ a_volts(1);
for k =2:length(a)
  c(k) =  a_volts(k-1)+a_volts(k);
end
c = c';
disp(c,'Duobinary coder output in volts:')
//Duobinary decoder output  by applying decision rule
for k =1:length(c)
  if(abs(c(k))>1)
    b_r(k) =  0;
  else
    b_r(k) = 1;
  end
end
b_r = b_r';
disp(b_r,'Recovered original sequence at detector oupupt:')
endfunction
//Result
//Precoder output in binary form:   
// 
//  1.    1.    0.    0.    1.    0.    0.  
// 
// Precoder output in volts:   
// 
//  1.    1.  - 1.  - 1.    1.  - 1.  - 1.  
// 
// Duobinary coder output in volts:   
//
//  2.    2.    0.  - 2.    0.    0.  - 2.  
// 
// Recovered original sequence at detector oupupt:   
// 
//  0.    0.    1.    0.    1.    1.    0.