DSPRelated.com
   Code Snippets

Computing Optimum Two-Stage Decimation Factors

Rick Lyons ● November 21, 2011 ● Coded in Matlab 
function [D1,D2] = Dec_OptimTwoStage_Factors(D_Total,Passband_width,Fstop)
%
%    When breaking a single large decimation factor (D_Total) into 
%    two stages of decimation (to minimize the orders of the 
%    necessary lowpass digital filtering), an optimum first 
%    stage of decimation factor (D1) can be found.  The second 
%    stage's decimation factor is then D_Total/D1.
%
%    Inputs:
%         D_Total = original single-stage decimation factor.
%         Passband_width = desired passband width of original 
%                          single-stage lowpass filter (Hz).
%         Fstop = desired beginning of stopband freq of original 
%                          single-stage lowpass filter (Hz).
%         
%    Outputs:
%         D1 = optimum first-stage decimation factor.
%         D2 = optimum second-stage decimation factor. 
%
%    Example: Assume you want to decimate a sequence by D_Total=100,
%             your original single-stage lowpass filter's passband 
%             width and stopband frequency are 1800 and 2200 Hz
%             respectively.  We use:
%             
%              [D1,D2] = Dec_OptimTwoStage_Factors(100, 1800, 2200)
%
%             giving us the desired D1 = 25, and D2 = 4. 
%             (That is, D1*D2 = 25*4 = 100 = D_Total.)
%
%    Author: Richard Lyons, November, 2011
	
% Start of processing
    DeltaF = (Fstop -Passband_width)/Fstop;
    
	Radical = (D_Total*DeltaF)/(2 -DeltaF); % Terms under sqrt sign.
	Numer = 2*D_Total*(1 -sqrt(Radical));  % Numerator.
	Denom = 2 -DeltaF*(D_Total + 1);  % Denominator.
	D1_estimated = Numer/Denom;  %Optimum D1 factor, but not 
                                  %                  an integer.

 % Find all factors of 'D_Total' 
	Factors_of_D_Total = Find_Factors_of_D_Total(D_Total);
 
 % Find the one factor of 'D_Total' that's closest to 'D1_estimated' 
    Temp = abs(Factors_of_D_Total -D1_estimated);
    Index_of_Min = find(Temp == min(Temp)); % Index of optim. D1
    D1 = Factors_of_D_Total(Index_of_Min); % Optimum first 
                                           %    decimation factor
    D2 = D_Total/D1;  % Second decimation factor.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allfactors] = Find_Factors_of_D_Total(X)
% Finds all integer factors of intger X.
% Filename was originally 'listfactors.m', written 
% by Jeff Miller.  Downloaded from Matlab Central.
	[b,m,n] = unique(factor(X));
%b is all prime factors
	occurences = [m(1) diff(m)];
	current_factors = [b(1).^(0:occurences(1))]';
	for index = 2:length(b)
        currentprime = b(index).^(0:occurences(index));
        current_factors = current_factors*currentprime;
        current_factors = current_factors(:);
	end
	allfactors = sort(current_factors);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Jesus Selva ● November 17, 2011 ● Coded in Matlab 
function Demo

while 1
  
  %Frequency grid with arbitrary start value, spacing, and number of elements. ...
  %Mathematically, the grid is
  %
  %   fo+Df*m, m=0, 1,...,M-1.
  %
  
  fo = 10*rand(1)-0.5; %Start value.
  Df = rand(1)+0.01;   %Frequency spacing. 
  M = round(rand(1)*1000); %Number of frequencies. 

  %Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
  %the grid is
  %
  %   to+Dt*n, n=0, 1,...,N-1.
  %

  to = 10*rand(1)-0.5; %Start value; (instant of first sample).
  Dt = rand(1)+0.01;   %Time spacing. 
  N = round(rand(1)*1000); %Number of samples. 
  
  %Random vector of samples.
  x = randn(1,N)+1i*randn(1,N);
  
  %x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N. 

  %We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.

  %Compute DFT using the direct and chirp-z methods.
  
  tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
  tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;
  
  disp(['Timing direct method: ', num2str(tmD), ' sec.'])
  disp(['Timing chirp z: ',num2str(tmF),' sec.'])
  disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
  disp(' ')
  input('(Press a key for another trial or Ctrl-C) ')
  disp(' ')
end

function X = chirpzDFT(x,to,Dt,fo,Df,M)

%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
%               N
%       X(k) = sum  x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
%              n=1
%
%Input parameters:
%
%  x: Signal samples. 
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
%  M: Desired number of output samples. 
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%

x = x(:).';

N = numel(x);
P = 2^ceil(log2(M+N-1));

%The next three lines do not depend on the vector x, and so they can be pre-computed if 
%the time and frequency grids do not change, (i.e, for computing the transform of 
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and 
%three if the grids change. 

a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse. 

%Weigh x using a and perform FFT convolution with phi. 
X = ifft( fft(x.*a,P) .* phi );

%Truncate the convolution tails and weigh using b. 
X = X(N:M+N-1) .* b;

function X = DirectGridDFT(x,to,Dt,fo,Df,M)

%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.

x = x(:).';

N = length(x);

X = zeros(1,M);

for k = 1:M
  for n = 1:N
    X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
  end
end

NCO bit true model in a single equation

kaz - ● October 23, 2011 ● Coded in Matlab 
%NCO bit true model
clear all; close all;

%example NCO parameters
n = 20000;                     %number of test samples
A = 2^15 - 1;                  %max amplitude, amplitude resolution 16 bit signed
M = 2^32;                      %NCO accumulator phase resolution
inc = round(M*.001);           %NCO accumulator phase increment 
k = 2^12;                      %lut phase resolution (lut size), may equal M

%single equation for NCO bit true model
y1 = round(A*sin(2*pi*round(k*mod((0:n-1)*inc/M,M))/k));      

%%%%%%%%%%%%%%%%%%%%%%%%% equivalent code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lsb = log2(M) - log2(k);                 %LSBs discarded when addressing
lut = round(A*sin(2*pi*(0:k-1)/k));      %lut, one cycle sine data

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

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

Complex filter with three multipliers per tap

kaz - ● October 22, 2011 ● Coded in Matlab 
%random input
I_in = randn(1,1024);
Q_in = randn(1,1024);

%random test filter
h = randn(1,20)+j*randn(1,20);    

%split up into 3 subfilters
h1 = real(h)+imag(h);         %extra adder
h2 = real(h);
h3 = real(h)-imag(h);         %extra subtractor

%apply subfilter
ch1 = filter(h1,1,I_in);
ch2 = filter(h2,1,Q_in-I_in); %extra subtractor
ch3 = filter(h3,1,Q_in);

%combine outputs
yr = round((ch3-ch2)*2^15);
yi = round((ch1+ch2)*2^15);
IQ_out = complex(yr,yi);

%direct computation
ref = round(2^15*filter(h,1,complex(I_in,Q_in)));
plot(IQ_out - ref);

Even/Odd FIR filter structure

kaz - ● October 22, 2011●1 comment ● Coded in Matlab 
% odd/even filter structure verification
% two filter structures compared
% (1) direct use of filter function
% (2) even/odd split structure

x = randn(1,2048);
h = fir1(23,.3);
y1 = filter(h,1,x); % direct filter for reference

he = h(1:2:end); % even coeff, relative to hardware definition
ho = h(2:2:end); % odd coeffs
xe = x(1:2:end); % even samples
xo = x(2:2:end); % odd samples

f1 = filter(ho,1,xe);   %subfilter 1
f2 = filter(he,1,xo);   %subfilter 2
f3 = filter(he,1,xe);   %subfilter 3
f4 = filter(ho,1,xo);   %subfilter 4

%add up subfilters, advancing f3 one stage
ye = [0 f1] + [0 f2];  
yo = [f3 0] + [0 f4];   

%interleave even/odd into one stream
y2 = [y1 0 0];
y2(1:2:end) = ye;
y2(2:2:end) = yo;

%quantise and compare both ref and subfilters
y1 = round(y1*2^15);
y2 = round(y2*2^15);
plot(y1 - y2(2:end-1))

Vector alignment for verification

kaz - ● October 22, 2011 ● Coded in Matlab 
%"align_run" is top level for testing "align" function

function align_run  
for i = 1:1
    L1 = 3400;   %length of first vector
    L2 = 4023;   %length of second vector
   
    %generate two random vectors
    x1 = randn(1,L1);
    x2 = randn(1,L2);
   
    %generate random index values for two identical segments
    r1 = round(rand(1)*L1);
    r2 = round(rand(1)*L2);
    r1(r1 == 0) = 1;
    r2(r2 == 0) = 1;
   
    %length of identical segments as percent of vector length
    n = round(50/100 * min([L1 L2]));
    n(n >= min([L1 L2])) = min([L1 L2])-1;
    while r1+n > L1, r1 = r1-1; end;
    while r2+n > L2, r2 = r2-1; end;
   
    %create identical segments at random location
    x1(r1:r1+n) = x2(r2:r2+n);
   
    %call align function
    [y1,y2] = align(x1,x2);
   
    %check alignment
    subplot(2,1,1);
    plot(y1);hold
    plot(y2,'r--') ;
    legend('y1','y2')
   
    subplot(2,1,2);
    plot(y1 - y2,'g')
    legend('(y1 - y2)')
   
end

return

%%%%%%%%%%%%%%%%% align %%%%%%%%%%%%%%%%%%
function [y1,y2] = align(x1,x2)

%force a column orientation
x1 = x1(:);
x2 = x2(:);

L1 = length(x1);
L2 = length(x2);
L = max([L1 L2]);

%equalise vector lengths
if L1 > L2, x2 = [x2; zeros(L1-L2,1)]; end;
if L1 < L2, x1 = [x1; zeros(L2-L1,1)]; end;

%apply correlation and find index of maximum output
y = xcorr(x1,x2);
n = find(y == max(y)) ;
n = mod(n,L)+1;

%align x1 vector with x2
y1 = [x1(n:end); x1(1:n-1)] ;
y2 = x2;

return

integrate RMS phase error from spectrum

Markus Nentwig ● October 21, 2011●1 comment ● Coded in Matlab 
% ****************************************************************
% RMS phase error from phase noise spectrum
% Markus Nentwig, 2011
%
% - integrates RMS phase error from a phase noise power spectrum
% - generates an example signal and determines the RMS-average
%   phase error (alternative method)
% ****************************************************************
function pn_snippet()
    close all;
    
    % ****************************************************************
    % Phase noise spectrum model
    % --------------------------
    % Notes: 
    % - Generally, a phase noise spectrum tends to follow
    %   |H(f)| = c0 + c1 f^n1 + c2 f^n2 + c3 f^n3 + ...
    %   Therefore, linear interpolation should be performed in dB
    %   on a LOGARITHMIC frequency axis.
    % - Single-/double sideband definition: 
    %   The phase noise model is defined for -inf <= f <= inf
    %   in other words, it contributes the given noise density at
    %   both positive AND negative frequencies. 
    %   Assuming a symmetrical PN spectrum, the model is evaluated 
    %   for |f| => no need to explicitly write out the negative frequency
    %   side.
    % ****************************************************************
    
    % PN model
    % first column: 
    %   frequency offset
    % second column: 
    %   spectral phase noise density, dB relative to carrier in a 1 Hz 
    %   observation bandwidth

    d = [0 -80
         10e3 -80
         1e6 -140
         9e9 -140];    
    f_Hz = d(:, 1) .';
    g_dBc1Hz = d(:, 2) .' -3;
    
    % get RMS phase error by integrating the power spectrum
    % (alternative 1)
    e1_degRMS = integratePN_degRMS(f_Hz, g_dBc1Hz)

    % get RMS phase error based on a test signal
    % (alternative 2)
    n = 2 ^ 20;
    deltaF_Hz = 2;    
    e2_degRMS = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz)
end

% ****************************************************************
% Integrate the RMS phase error from the power spectrum
% ****************************************************************
function r = integratePN_degRMS(f1_Hz, g1_dBc1Hz)
    1;
    % integration step size
    deltaF_Hz = 100;
    
    % list of integration frequencies
    f_Hz = deltaF_Hz:deltaF_Hz:5e6;
    
    % interpolate spectrum on logarithmic frequency, dB scale
    % unit is dBc in 1 Hz, relative to a unity carrier
    fr_dB = interp1(log(f1_Hz+eps), g1_dBc1Hz, log(f_Hz+eps), 'linear');
    
    % convert to power in 1 Hz, relative to a unity carrier
    fr_pwr = 10 .^ (fr_dB/10);

    % scale to integration step size
    % unit: power in deltaF_Hz, relative to unity carrier
    fr_pwr = fr_pwr * deltaF_Hz;

    % evaluation frequencies are positive only
    % phase noise is two-sided 
    % (note the IEEE definition: one-half the double sideband PSD)    
    fr_pwr = fr_pwr * 2;

    % sum up relative power over all frequencies
    pow_relToCarrier = sum(fr_pwr);

    % convert the phase noise power to an RMS magnitude, relative to the carrier
    pnMagnitude = sqrt(pow_relToCarrier);

    % convert from radians to degrees
    r = pnMagnitude * 180 / pi;
end

% ****************************************************************
% Generate a PN signal with the given power spectrum and determine
% the RMS phase error
% ****************************************************************
function r = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz);
    A = 1; % unity amplitude, arbitrary

    % indices of positive-frequency FFT bins. 
    % Does not include DC
    ixPos = 2:(n/2);
    
    % indices of negative-frequency FFT bins. 
    % Does not include the bin on the Nyquist limit
    assert(mod(n, 2) == 0, 'n must be even');
    ixNeg = n - ixPos + 2;
    
    % Generate signal in the frequency domain
    sig = zeros(1, n);
    % positive frequencies
    sig(ixPos) = randn(size(ixPos)) + 1i * randn(size(ixPos));
    % for purely real-valued signal: conj()
    % for purely imaginary-valued signal such as phase noise: -conj()
    % Note: 
    %   Small-signals are assumed. For higher "modulation indices", 
    %   phase noise would become a circle in the complex plane
    sig(ixNeg) = -conj(sig(ixPos));
    
    % normalize energy to unity amplitude A per bin
    sig = sig * A / RMS(sig);

    % normalize energy to 0 dBc in 1 Hz
    sig = sig * sqrt(deltaF_Hz);

    % frequency vector corresponding to the FFT bins (0, positive, negative)
    fb_Hz = FFT_freqbase(n, deltaF_Hz);

    % interpolate phase noise spectrum on frequency vector
    % interpolate dB on logarithmic frequency axis
    H_dB = interp1(log(f_Hz+eps), g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');

    % dB => magnitude
    H = 10 .^ (H_dB / 20);

    % plot on linear f axis
    figure(); hold on;
    plot(fftshift(fb_Hz), fftshift(mag2dB(H)), 'k');
    xlabel('f/Hz'); ylabel('dBc in 1 Hz');
    title('phase noise spectrum (linear frequency axis)');
    
    % plot on logarithmic f axis
    figure(); hold on;
    ix = find(fb_Hz > 0);
    semilogx(fb_Hz(ix), mag2dB(H(ix)), 'k');
    xlabel('f/Hz'); ylabel('dBc in 1 Hz');
    title('phase noise spectrum (logarithmic frequency axis)');
    
    % apply frequency response to signal
    sig = sig .* H;
    
    % add carrier
    sig (1) = A;

    % convert to time domain
    td = ifft(sig);
    
    % determine phase 
    % for an ideal carrier, it would be zero
    % any deviation from zero is phase error
    ph_deg = angle(td) * 180 / pi;

    figure(); 
    ix = 1:1000;
    plot(ix / n / deltaF_Hz, ph_deg(ix));
    xlabel('time / s');
    ylabel('phase error / degrees');
    title('phase of example signal (first 1000 samples)');
    
    figure(); 
    hist(ph_deg, 300);
    title('histogram of example signal phase');
    xlabel('phase error / degrees');

    % RMS average
    % note: exp(1i*x) ~ x
    %       as with magnitude, power/energy is proportional to phase error squared
    r = RMS(ph_deg);
    
    % knowing that the signal does not have an amplitude component, we can alternatively
    % determine the RMS phase error from the power of the phase noise component

    % exclude carrier:
    sig(1) = 0; 

    % 3rd alternative to calculate RMS phase error
    rAlt_deg = sqrt(sig * sig') * 180 / pi    
end

% gets a vector of frequencies corresponding to n FFT bins, when the 
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
    fb_Hz = 0:(n - 1);
    fb_Hz = fb_Hz + floor(n / 2);
    fb_Hz = mod(fb_Hz, n);
    fb_Hz = fb_Hz - floor(n / 2);
    fb_Hz = fb_Hz * deltaF_Hz;
end

% Root-mean-square average
function r = RMS(vec)
    r = sqrt(vec * vec' / numel(vec));
end

% convert a magnitude to dB
function r = mag2dB(vec)
    r = 20*log10(abs(vec) + eps);
end

Saving a Matlab figure to file using command-line

Jeff T ● September 30, 2011 ● Coded in Matlab 
function save_fig(h, name, format)
% Matlab has a wierd way of saving figures and keeping the plot looking
% exactly the way it does on your screen.  Thi function simplifies the 
% entire process in a way that is barely documented. 
%
% Usage:      SAVE_FIG(h, name, format);
%
%             H is the handle to the figure.  This can be obtain in the 
%               following manner:  H = figure(1);
%             NAME is the filename without extension
%             FORMAT is the graphic format.  Options are:
%             
%                   'bmpmono'    BMP monochrome BMP
%                   'bmp16m'     BMP 24-bit BMP
%                   'bmp256'     BMP 8-bit (256-color) 
%                   'bmp'        BMP 24-bit	
%                   'meta'       EMF
%                   'eps'        EPS black and white
%                   'epsc'       EPS color
%                   'eps2'       EPS Level 2 black and white
%                   'epsc2'      EPS Level 2 color
%                   'fig'        FIG Matlab figure file
%                   'hdf'        HDF 24-bit
%                   'ill'        ILL (Adobe Illustrator)
%                   'jpeg'       JPEG 24-bit
%                   'pbm'        PBM (plain format) 1-bit
%                   'pbmraw'     PBM (raw format) 1-bit
%                   'pcxmono'    PCX 1-bit
%                   'pcx24b'     PCX 24-bit color PCX file format
%                   'pcx256'     PCX 8-bit newer color (256-color)
%                   'pcx16'      PCX 4-bit older color (16-color)
%                   'pdf'        PDF Color PDF file format
%                   'pgm'        PGM Portable Graymap (plain format)
%                   'pgmraw'     PGM Portable Graymap (raw format)
%                   'png'        PNG 24-bit
%                   'ppm'        PPM Portable Pixmap (plain format)
%                   'ppmraw'     PPM Portable Pixmap (raw format)
%                   'svg'        SVG Scalable Vector Graphics
%                   'tiff'       TIFF 24-bit
%
% Author: sparafucile17

if(strcmp(format,'fig'))
    saveas(h, name, 'fig');
else
    options.Format = format;
    hgexport(h, name, options);
end

Putting a Matlab figure to a specific screen location

Jeff T ● September 30, 2011 ● Coded in Matlab 
function set_fig_position(h, top, left, height, width)
% Matlab has a wierd way of positioning figures so this function
% simplifies the poisitioning scheme in a more conventional way.
%
% Usage:      SET_FIG_POISITION(h, top, left, height, width);
%
%             H is the handle to the figure.  This can be obtain in the 
%               following manner:  H = figure(1);
%             TOP is the "y" screen coordinate for the top of the figure
%             LEFT is the "x" screen coordinate for the left side of the figure
%             HEIGHT is how tall you want the figure
%             WIDTH is how wide you want the figure
%
% Author: sparafucile17

% PC's active screen size
screen_size = get(0,'ScreenSize');
pc_width  = screen_size(3);
pc_height = screen_size(4);

%Matlab also does not consider the height of the figure's toolbar...
%Or the width of the border... they only care about the contents!
toolbar_height = 77;
window_border  = 5;

% The Format of Matlab is this:
%[left, bottom, width, height]
m_left   = left + window_border;
m_bottom = pc_height - height - top - toolbar_height - 1;
m_height = height;
m_width  = width - 1;

%Set the correct position of the figure
set(h, 'Position', [m_left, m_bottom, m_width, m_height]);

%If you want it to print to file correctly, this must also be set
% and you must use the "-r72" scaling to get the proper format
set(h, 'PaperUnits', 'points');
set(h, 'PaperSize', [width, height]); 
set(h, 'PaperPosition', [0, 0, width, height]); %[ left, bottom, width, height]

Least Squares fit to a linear equation

Jeff T ● September 30, 2011 ● Coded in Matlab 
function [m, b] = ls_linear(x_observed, y_observed)
%
% Perform Least Squares curve-fitting techniques to solve the coefficients
% of the linear equation: y = m*x + b.  Input and Output equation 
% observations must be fed into the algorithm and a best fit equation will
% be calculated.
%
% Usage:     [M,B] = ls_linear(observed_input, observed_output);
%
%            observed_input  is the x-vector observed data
%            observed_output is the y-vector observed data
%
% Author:    sparafucile17 06/04/03
%

if(length(x_observed) ~= length(y_observed))
    error('ERROR: Both X and Y vectors must be the same size!');
end

%Calculate vector length once
vector_length = length(x_observed);

%
% Theta = [   x1    1   ]
%         [   x2    1   ]
%         [   x3    1   ]
%         [   x4    1   ]
%         [   ...  ...  ]
%
Theta = ones(vector_length, 2);
Theta(1:end, 1) = x_observed;

%
% Theta_0 = [  y1   ]
%           [  y2   ]
%           [  y3   ]
%           [  y4   ]
%           [  ...  ]
%
Theta_0 = y_observed;

%
% Theta_LS = [  M  ]
%            [  B  ]
%
Theta_LS = ((Theta' * Theta)^-1 ) * Theta' * Theta_0;

%Return the answer
m = Theta_LS(1);
b = Theta_LS(2);