DSPRelated.com
   Code Snippets

passive Sum difference method for finding bearing

● December 15, 2010 ● Coded in Matlab 
clc
clear all
close all

fs=20e7;
T=1/fs;

t=0:T:0.0000002;
f=10e6;
scanagle_deg=60; % Range covered by seeker antenna
step_deg=6;

% for ploting original arrays of sum, diff and ratio
theta=-deg2rad(scanagle_deg):deg2rad(step_deg):deg2rad(scanagle_deg);

for i=1:length(theta);

ant1=1*sin(2*pi*f*t);               %Antenna 1
ant2=1*sin(2*pi*f*t+theta(i));      %Antenna 2

[my,mx]=max(ant1);

sum_ant(i)=ant1(mx)+ant2(mx);               %Sum of antennas
diff_ant(i)=ant1(mx)-ant2(mx);              %diff of antennas
ratio_ant(i)=diff_ant(i)/sum_ant(i);        %ratio

end

% subplot(311)
% plot(t,ant1,t,ant2)
% subplot(211)
% plot(rad2deg(theta),sum_ant,rad2deg(theta),diff_ant)
% subplot(212)
% plot(rad2deg(theta),ratio_ant)

% 
% figure
% subplot(211)
% plot(rad2deg(theta),sumn_ant,rad2deg(theta),diffn_ant)
% subplot(212)
% plot(rad2deg(theta),ration_ant)

%%%%%%%%%%%%%%%%Recieved signal at antenna for DOA%%%%%%%%%%%%%%%%%%%%%%%%%%

A_ant1=2;                   % Amplitude of antenna 1
A_ant2=3;                   % Amplitude of antenna 2

DOA_angle_target=20;        % DOA of TARGET

sim_ant11=A_ant1*sin(2*pi*f*t);                         % Simulated antenna 1
sim_ant22=A_ant2*sin(2*pi*f*t+deg2rad(DOA_angle_target));   % Simulated antenna 2

sim_ant1=sim_ant11/max(sim_ant11);          %normalization
sim_ant2=sim_ant22/max(sim_ant22);

[smy,smx]=max(sim_ant1);

%%%%also take for same sample value of data

sim_sum_ant=sim_ant1(smx)+sim_ant2(smx);
sim_diff_ant=sim_ant1(smx)-sim_ant2(smx);
sim_ratio_ant=sim_diff_ant/sim_sum_ant;

%%%%%%%%%%%%% Polynomial fitting of ratio_ant of 6th order%%%%%%%%%%%%
% Error in DOA is obtained because of this
% if polynomial fitting is removed or more accurate results are obtained
% DOA of target can be found with greater accuracy
% Polynomial was obtained though curve fitting toolbox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p1 =  1.047e-008;
p2 = -6.907e-007;
p3 =  2.383e-005;
p4 =  -0.0004914;
p5 =    0.008555;
p6 =    -0.09626;
p7 =      0.4215;
x=0;

for ii=1:2200
    x=x+0.01;
fitted_model_rat_ant(ii) = p1*x^6 + p2*x^5 + p3*x^4 + p4*x^3 + p5*x^2 + p6*x + p7;

end

theta_new=-scanagle_deg:.0545:scanagle_deg;    % Theta for ploting fitted model
figure
plot(theta_new(1:2200),fitted_model_rat_ant)

c1=ceil(sim_ratio_ant*7000);           % For comparison of sim_ratio ant and model
c2=ceil(fitted_model_rat_ant*7000);    % Different threshold Threshold can be chosen 
[r,c,v]= find(c2==c1);

detected_theta=(c.*0.0545)-60       %theta from curve fitting model

if(A_ant1>A_ant2)           % condition for checking which angle was correct
    correct_theta=detected_theta(1)
elseif(A_ant1<A_ant2)
    correct_theta=detected_theta(2)
else
     correct_theta=detected_theta
end

Goertzel Filterbank to the Implementation of a Nonuniform DFT

Miguel De Jesus Rosario ● December 14, 2010●1 comment ● Coded in Matlab 
% VectorGoertzel    Goertzel's Algorithm filter bank.
%   Realization of the Goertzel's Algorithm to compute the Nonuniform DFT 
%   of a signal(a column vector named signalw) of length Nw with sampling 
%   frecuency fs at the desired frecuencies contained in vector f. The 
%   function returns the NDFT magnitude in a vector with the same length of f.

function xk=Gfilterbank(signalw,f,fs,Nw)

% Inititialization of the different variables 
n=2;
signalw=signalw(1:Nw);

cost=cos(2*pi*f/fs);
sint=sin(2*pi*f/fs);
L=length(cost);

y1(1:L)=0;
y2(1:L)=0;
signalw=[0 0 signalw]; %Signal is delayed by two samples

% Goertzel Feedback Algorithm

while((n-2) < Nw)
    n=n+1;
    xnew(1:L)=signalw(n);
    y=xnew+2*cost.*y1-y2;
    y2=y1;
    y1=y;
end

% Goertzel Forward Algorithm

rey=y1-y2.*cost;
imy=y2.*sint;

% Magnitude Calculation

xk=abs(rey+imy*j)';

Discrete Wavelet Transform - Chain processing (Linear Convolution)

David Valencia ● December 5, 2010●1 comment ● Coded in Matlab 
% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% Author: David Valencia
% UPIITA IPN 2010
%
% Posted in DSPrelated.com
% http://www.dsprelated.com/showcode/47.php
%
% Computes the Discrete Wavelet Transform of
% n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% chain processing method (linear convolution)
%
% Dependencies: 
%   formfilters.m - http://www.dsprelated.com/showcode/44.php
%   formfiltersdwt.m - http://www.dsprelated.com/showcode/45.php
%   upsample2.m - http://www.dsprelated.com/showcode/10.php
%   recorreder.m - http://www.dsprelated.com/showcode/43.php
%   getsincdwt.m - http://www.dsprelated.com/showcode/46.php
%
% Revisions:
%       v1.0a Commented and translated in English
% ----------------------------------------------------

close all; clear all; clc;

% ######################################

disp('== CHAIN PROCESSING DWT ==')

% Filter parameters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db2';

if(typeofbasis == 'b')
    [Rf,Df] = biorwavf(typbior);
    [h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
    [h0,h1,g0,g1] = wfilters(typor);
end;

% Input values, change as needed
N = 16;
x = (1:N);
L = length(h0);

figure(1);
subplot(2,1,1);
stem(x);
xlabel('Original signal');

% Filter bank parameters, change as needed
n_stages = 3; %Of the low-pass part
% Checked for the minimal case: DWT = DWPT n_stages = 1
n_branches = n_stages + 1; %Of the low-pass part
niter = 300; %Number of iterations

% Define input buffer dimensions
hx = formfilters(n_stages, 1, h0, h1);
lhx = length(hx);

% Input buffer vector
inputbuffer = zeros(lhx,1);

% Synthesis input buffer
sbuffer = zeros(n_branches,lhx);

% This vector will store the length of each branch filter
Lfiltros = zeros(n_branches, 1);

% This vector will store the decimate factor for each branch
dec_factors = zeros(n_branches, 1);

% This vector will store the synchronization factor for each branch
sinc_factors = zeros(n_branches, 1);

% Cycle to get the synchronization factors
for j = 0:n_branches-1
    try
        sinc_factors(j+1) = getsincdwt(n_stages, j+1, L);
    catch error
        disp('ERROR: Error, experimentation is needed for this case');
        disp('The results may not be correct');
        % return;
    end
end

%% Filter matrices formation
% This matrix will store each of the branch filters as vectors
H = zeros(n_branches, lhx); %Malloc
G = zeros(n_branches, lhx); %Malloc

for i = 0:n_stages
    
    if i >= n_stages-1
        dec_factors(i+1) = 2^n_stages;
    else
        dec_factors(i+1) = 2^(i+1);
    end
    
    hx = fliplr(formfiltersdwt(n_stages,n_stages+1-i,h0,h1));
    gx = fliplr(formfiltersdwt(n_stages,n_stages+1-i,g0,g1));
    lhx = length(hx);
    lgx = length(gx);
    Lfiltros(i+1) = lhx;
    H(i+1,1:lhx) = hx;
    G(i+1,1:lhx) = gx;
end

% Debug code
disp('Dimensiones de los filtros')
Lfiltros
disp('Factores de diezmado')
dec_factors
disp('Matriz de filtros análisis: ');
H
disp('Matriz de filtros síntesis: ');
G

chanbufftestigo = zeros(n_branches,1);
outputbuffer = zeros(1,niter);

%% MAIN CYCLE
for i = 1:niter %General FOR (for iterations)
    i % Print iteration number
    % Shifting of input buffer (sequentially)
    
    inputbuffer = recorreder(inputbuffer,1);
    if i<=N
        inputbuffer(1)=x(i);
    else
        inputbuffer(1)=0;
    end
    
    inputbuffer %Print buffer (debug)
    
    %% Analyisis stage (h0 and h1)
    % The following cycle will select each of the branches for processing
    % If the iteration counter is a multiply of the decimation factor,
    % the convolution is calculated, otherwise, zeros are sent to the
    % channel
    
    y = zeros(n_branches,1); %Clear the output buffer
    
    for j = 0:n_branches-1
        fprintf('\nBranch: %d, Decimating factor: %d\n',j+1,dec_factors(j+1));
        fprintf('\nFilter length: %d\n',Lfiltros(j+1));
        if mod(i,dec_factors(j+1)) == 1
            fprintf('i = %d is a multiply of the factor: %d\n',i,dec_factors(j+1));
            tempfilt = H(j+1,1:Lfiltros(j+1));
            % If the current iteration (clock cycle) is a multiply of the
            % factor, the convolution is computed
            for k=1:Lfiltros(j+1) %convolution
                y(j+1,1) = y(j+1,1) + inputbuffer(k)*tempfilt(k);
            end;
        end
        disp('Results in the channel');
        chanbuff(:,1) = y %Debug printing
    end
    
    %% Synthesis stage (g0 and g1)
    
    % Shifting and interpolation of the synthesis buffer
    for j = 1:n_branches
        sbuffer(j,:) = recorreder(sbuffer(j,:),1);
        % Interpolation of the synthesis stage inputs
        if mod(i,dec_factors(j)) == 1
            
            fprintf('Inserting sample to the synthesis branch %d with dec_factor %d\n',j,dec_factors(j));
            sbuffer(j,1) = chanbuff(j,1);
        else
            fprintf('Inserting a zero to the synthesis branch %d with dec_factor %d\n',j,dec_factors(j));
            sbuffer(j,1) = 0;
        end
    end
    
    disp('Synthesis buffer matrix');
    sbuffer
    
    xnbuff = zeros(n_branches,1);
    
    for j=0:n_branches-1
        fprintf('Branch: %d , dec_factor = %d',j+1,dec_factors(j+1));
        
        % === BRANCH SYNCHRONIZATION ===
        % The branches (except the two last ones down the bank filter)
        % will only take a part of the vector that corresponds with their
        % channel, taking 'Lx' elements, being Lx the number of coefficient
        % that their filter has. An offset is applied
        
        if j < n_branches - 2
            [m,n] = size(sbuffer);
            tempsinth = sbuffer(j+1,n-Lfiltros(j+1)+1:n) %TESTING
        else % The lowest branches take the whole vector
            tempsinth = sbuffer(j+1,1+sinc_factors(j+1):Lfiltros(j+1)+sinc_factors(j+1))
        end
        
        for k = 1 : Lfiltros(j+1) %Convolution
            xnbuff(j+1,1) = xnbuff(j+1,1) + tempsinth(k)*G(j+1,k);
        end
    end
    
    xnbuff
    outputbuffer(i) = sum(xnbuff);
    disp('Output buffer: ');
    if i > N
        disp(outputbuffer(1,i-N+1:i))
        figure(1);
        subplot(2,1,2);
        stem(outputbuffer(1,i-N+1:i));        
        xlabel('Recovered signal');
    else
        disp(outputbuffer(1,1:N))
        figure(1);
        subplot(2,1,2);
        stem(outputbuffer(1,1:N));
        xlabel('Recovered signal');        
    end
    pause; 
end

% END OF PROGRAM >>>

getsincdwt.m - Get DWT synchronization for chain-processing

David Valencia ● December 5, 2010 ● Coded in Matlab 
% UPIITA IPN 2010
% Valencia Pesqueira José David

% Function used to obtain the synchronization factor for each
% branch in the DWT chain-processing program

function [x] = getsincdwt(n_etapas, rama, Lfiltrobase)

switch(n_etapas)
    case 1 %One stage, two-branches. DWPT equivalent
        x = 0;
        
    case 2  %2 stages, 3 branches
        switch(rama)
            case 1 % High-pass branch
                x = (Lfiltrobase-1)*2;
            otherwise
                x = 0;
        end
        
    case 3 %3 branches, 4 stages
        switch(rama)
            case 1 % High-pass branch
                x = (Lfiltrobase-1)*2;
            case 2
                x = 10;
            otherwise
                x = 0;
        end
end

% Experimental ANALISIS
% 2 etapas
% lfiltro - sinc_factor experimental
% 
% 4 - 6       (4-1)*2 = 6
% 6 - 10	    (6-1)*2 = 10
% 8 - 14	    (8-1)*2 = 7

formfiltersdwt.m DWT filter generator adaptor

David Valencia ● December 5, 2010 ● Coded in Matlab 
% By David Valencia
% As posted in DSPRelated.com
% at http://www.dsprelated.com/showcode/45.php
% Used for the DWT chain processing program

function [hx] = formfiltersdwt(n_stages,branch,h0,h1)
if(branch==1)
    hx=formfilters(n_stages,branch,h0,h1);
elseif(branch==2)
    hx=formfilters(n_stages,2,h0,h1);
else
    hx=formfilters(n_stages-branch+2,2,h0,h1);
end

% The code line that is changed in formfilters.m for DWT is:
%p = mod(2^(n_stages-1-i),p);

DWT filter generator - formfilters.m (for chain processing)

David Valencia ● December 5, 2010 ● Coded in Matlab 
% Created by José David Valencia Pesqueira
% UPIITA-IPN 2010
% For the DWT chain-processing example
% as posted in DSPrelated.com
% http://www.dsprelated.com/showcode/44.php
%
function [hx] = formfilters(n_stages,branch,h0,h1)
p = branch;

% Seed vector
hx = 0;
hx(1) = 1;

switch n_stages
    case 1
        if mod(branch,2) ~= 0
            hx = h0;
        else
            hx = h1;
        end
    case 2
        switch branch
            case 1
                hx = conv(h0,upsample2(h0,2));
            case 2
                hx = conv(h0,upsample2(h1,2));
            case 3
                hx = conv(h1,upsample2(h0,2));
            case 4
                hx = conv(h1,upsample2(h1,2));
            otherwise
                beep;
                fprintf('\nFor a 2-stage bank there cannot be a fifth branch');
        end
        
    otherwise
        
        for i=0:n_stages-2
            q = floor(p /(2^(n_stages-1-i)));
            if (q == 1)
                hx = conv(hx,upsample2(h1,2^i));
            else
                hx = conv(hx,upsample2(h0,2^i));
            end
            % p = mod(p,2^(n_stages-1-i)); %For DWPT
            p = mod(2^(n_stages-1-i),p); %For DWT
        end
        
        t = mod(branch,2);
        if(t == 1)
            hx = conv(hx,upsample2(h0,2^(n_stages-1)));
        else
            hx = conv(hx,upsample2(h1,2^(n_stages-1)));
        end
             
end

recorreder.m - Vector shifting made simple

David Valencia ● December 5, 2010 ● Coded in Matlab 
function [xd] = recorreder(x,N)
% This function shifts a vector by N positions to the right
% As posted in DSPrelated.com
% http://www.dsprelated.com/showcode/43.php
%
lx = length(x);
if N>lx
    fprintf('N has to be smaller than the vector length');
    return;
else
    xd(N+1:lx) = x(1:lx-N);
end;

Echo Filter

Miguel De Jesus Rosario ● December 1, 2010●3 comments ● Coded in Matlab 
% function echo_filter:
%     1) computes and returns the impulse response h of the LTI system (echo filter)
%     2) Filter the input signal.  The output signal is yecho. 
%
% INPUTS:
% signal = input signal
% times = vector containing the time delays (in seconds) of the repetitions(echoes)
% attenuations = vector containing the respective attenuations of the echoes
% fs = sampling frequency of the input signal.
%
% OUTPUTS:
% yecho = output signal (passed through the echo filter)
% h = impulse response of the echo filter.

function [yecho, h] = echo_filter(signal, times, attenuations, fs)

h(1)=1; % This is the first coefficient of the echo filter

% Computing the impulse response h
for i=1:length(times),
    samples(i) = times(i)*fs; % Calculating the sample-times (N = t*fs)
    h(floor(samples(i))) = attenuations(i); % Impulse response coeficients
end

% #########################################################################
% You may use the following implementation instead of the illustrated
% above:

%       samples = times*fs;
%       h(floor(samples)) = attenuations;

% In this case, the implementation is done without a for-loop.
% #########################################################################

yecho = filter(h,1,signal(:,1)); % filtering of the input signal results in a signal with echoes

% You may test the filter using the Matlab signal "splat" (write load splat
% in the Matlab command window).  Use function sound() in order to
% listening both, the input and the output, signals.

Beampattern of a Linear Array of Antennas (Array Processing)

Miguel De Jesus Rosario ● November 30, 2010●1 comment ● Coded in Matlab 
% The function NDFT2 computes the sinc pattern of a linear array of 
% sensors. The function receives three input: the interelement spacing d,
% the array weigths a, and the zero padding desired Npad.

function NDFT2(d, a, Npad)

k  = -Npad/2 : Npad/2-1; % index
u  = pi*k/Npad; % u is a vector with elements in the range of -pi/2 to pi/2
uk = sin(u);
n  = 0:Npad-1; % n is an index
m  = 1:Npad; % m is an index

% Wavenumber K = 2*pi/landa 
% d = is a fraction or a multiple of lambda.
% v = K*d*uk(m).
  
v = 2*pi*d*uk(m)';
Dn(m,n+1) = exp(j*v*n);

puk = Dn * a'; % Computing the Beampattern

% Plot of the Beampatterns
figure(1); subplot(2,1,1); plot(uk,20*log10(abs(puk)));grid on; xlabel('sin(theta)'); ylabel('Magnitude (dB)')
axis([-1 1 -100 0])
subplot(2,1,2); plot(uk,puk);grid on; xlabel('sin(theta)'); ylabel('Magnitude')
axis([-1 1 -1 1])
warning off;
% Plot the beampattern as a polar graph
figure(2); polar(u',puk); hold on; polar(-u',puk);hold off         

%  Function pattern()
%
%     The function pattern() computes and plots the beampattern of a 
%     Linear Array of sensors. This function has three inputs: the number of elements in
%     the array, the pointing direction of the beam pattern (an angular
%     direction in radians), and a constant spacing between the elements of 
%     the array (as fraction of the wavelenght(lambda)of the transmitted  
%     signal.  The optimum spacing between elements of the array is 0.5. 
%     You must also select any of the windows presented in the menu. 
%     Windowing techniques are used to reduced the sidelobes of the pattern.   
%     The list of available windows is the following:
%
%           @bartlett       - Bartlett window.
%           @barthannwin    - Modified Bartlett-Hanning window. 
%           @blackman       - Blackman window.
%           @blackmanharris - Minimum 4-term Blackman-Harris window.
%           @bohmanwin      - Bohman window.
%           @chebwin        - Chebyshev window.
%           @flattopwin     - Flat Top window.
%           @gausswin       - Gaussian window.
%           @hamming        - Hamming window.
%           @hann           - Hann window.
%           @kaiser         - Kaiser window.
%           @nuttallwin     - Nuttall defined minimum 4-term Blackman-Harris window.
%           @parzenwin      - Parzen (de la Valle-Poussin) window.
%           @rectwin        - Rectangular window.
%           @tukeywin       - Tukey window.
%           @triang         - Triangular window.
%
%     For example:      pattern(21, pi/4, 0.5);
%           Plots the beampattern of an linear array with 21 elements equally spaced  
%           at a half of the wavelenght(lambda/2), and a pointing
%           direction of 45 degrees. For uniform arrays use a rectangular 
%           window (rectwin).

function pattern(array_number, angular_direction, spacing)

close all
clc

N = array_number;
delta = angular_direction;
d = spacing;

Npad=1024;
n=0:N-1;

delta = 2*pi*d*sin(delta);
an=1/N*exp(j*n*delta);

for i=0:500000,
    
option = menu('Choose the desired Window', 'Bartlett', 'Barthannwin', 'Blackman', 'Blackmanharris', 'Bohmanwin', 'Chebwin', 'Flattopwin', 'Gausswin', 'Hamming', 'Hann', 'Kaiser', 'Nuttallwin', 'Parzenwin', 'Rectwin', 'Tukeywin', 'Triang', 'Exit'); 

switch option
    
    case 1
        close all
        clear a;
        a=an;
        a = a.*bartlett(N)';
        a(Npad)=0;
      
        NDFT2(d, a, Npad);
        
        
    case 2
        close all
        clear a;
        a=an;
        a = a.*barthannwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
        
    case 3
        close all 
        clear a;
        a=an;
        a = a.*blackman(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 4
        close all  
        clear a;
        a=an;
        a = a.*blackmanharris(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 5
        close all 
        clear a;
        a=an;
        a = a.*bohmanwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 6
        close all 
        clear a;
        a=an;
        a = a.*chebwin(N,40)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 7
        close all
        clear a;
        a=an;
        a = a.*flattopwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 8
        close all 
        clear a;
        a=an;
        a = a.*gausswin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 9
        close all  
        clear a;
        a=an;
        a = a.*hamming(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 10
        close all 
        clear a;
        a=an;
        a = a.*hann(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 11
         close all
         clear a;
        a=an;
        a = a.*kaiser(N,1)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 12
         close all
         clear a;
        a=an;
        a = a.*nuttallwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 13
         close all 
         clear a;
        a=an;
        a = a.*parzenwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 14
         close all 
         clear a;
        a=an;
        a = a.*rectwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 15
         close all  
         clear a;
        a=an;
        a = a.*tukeywin(N,0)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 16
        close all
        clear a;
        a=an;
        a = a.*triang(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 17
        
        break;
              
    end
    
end

WCDMA channelization code generator

Markus Nentwig ● November 13, 2010●1 comment ● Coded in Matlab 
% WCDMA channelization codes 
% source:
% 3GPP TS 25.213 V10.0.0 section 4.3.1.1
%
% parameters: 
% sf: spreading factor
% k: code number
%
% The returned code is a column vector with length sf
%
% this code has not been tested extensively. Please verify
% independently that it does what it promises.
function code=UTRA_FDD_chanCode(sf, k)
persistent flag;
persistent codes;

% * ********************************************
% * Build codebook
% * ********************************************
if isempty(flag)
  codes={};
  flag=1;
  
  % start of recursion: Identity code for sf=1
  codes{1, 1}=1;
  for sfi=1:8
    sfg=2 ^ sfi;
    for kgDest=0:2:sfg-2
      kgSrc=kgDest/2;
      prev=codes{sfg/2, kgSrc+1};
      % generation method:
      % - duplicate a lower-sf code
      % - duplicate and change sign
      codes{sfg, kgDest+1}=[prev prev];
      codes{sfg, kgDest+2}=[prev -prev];
    end
  end
end

% * ********************************************
% * look up the requested code from codebook
% * ********************************************
code=transpose(codes{sf, k+1});

% ################## CUT HERE ######################### 

% Example use (put this into a separate file)
sf=128; codenum=3;
chanCode=UTRA_FDD_chanCode(sf, codenum);

sig=[1 0 0 1 0 0 ]; % signal, row vector
len=size(sig, 2), 
% convolve:
s=chanCode * sig; % now matrix, one row per repetition
len=len*sf;
s=reshape(s, 1, len);
% plot
stem(s);