Code covered by the BSD License  

Highlights from
Bessel Polynomial Linearization Coefficients

image thumbnail
from Bessel Polynomial Linearization Coefficients by Robert Dickson
Linearization coefficients for Bessel polynomials, and certain rational polynomial spectra.

bpqpsdg(gi,f) 
function [go,psd] = bpqpsdg(gi,f)
%BPQPSDG  Power spectral density of a Gamma-Gamma process.
%
%   [go,psd] = bpqpsdg(gi,f)   Power spectral density psd of Gamma-Gamma
%   random process G with parameters specified by structure gi. PSD is
%   evaluated at vector f of frequencies. Empty elements of gi are
%   completed in go. 
%
%   gi = struct ( ...
%     'Fscale', 1, ...  % Frequency scale 
%     'Ratio', 1/8, ... % Ratio aX/aY
%     'alpha', 5, ...   % Gamma shape parameter of X
%     'beta', 1.5, ...  % Gamma shape parameter of Y
%     'nX', 4, ...      % K-Bessel spectral order underlying X
%     'nY', 4, ...      % K-Bessel spectral order underlying Y
%     'aX', [], ...     % Spectral scale factor of X
%     'aY', [], ...     % Spectral scale factor of Y
%     'SI', [], ...     % Scintillation Index
%     'sm2', [], ...    % Second spectral moment of G
%     'sm4', [] ...     % Fourth spectral moment of G
%   );
%
%   Notes:
%       SI =  1/alpha + 1/beta + 1/(alpha*beta)
%       sm2 = (2*pi*Fscale)^2
%
%   References:
%   [1] M.A. Al-Habash, L.C. Andrews, and R.L. Phillips, Mathematical model
%   for the irradiance probability density function of a laser beam
%   propagating through turbulent media, Optical Engineering,
%   Vol. 40 No. 8, August 2001.
%
%   Uses: bpqpsdx, bpqpsdw.
%

% 2013-06-11 Robert Dickson

if rem(gi.nX,1) ~= 0 || rem(gi.nY,1) ~= 0
    error('bpqpsdg:arg','nX,nY must be integer-valued scalars.');
end

minn = min(gi.nX,gi.nY);

if minn < 2
    error('bpqpsdg:arg','nX and nY must greater than or equal to two.');
end

go = gi;
go.sm4 = Inf;

go.SI = 1/go.alpha+1/go.beta+1/(go.alpha*go.beta); % Scintillation Index
CX = 1/(go.alpha*go.SI);
CY = 1/(go.beta*go.SI);
CXY = 1/(go.alpha*go.beta*go.SI);

go.sm2 = (2*pi*go.Fscale)^2;  % Second Spectral Moment
kr2 = 2*(2*go.nX-3)*(2*go.nY-3)*(1+go.alpha+go.beta)/ ...
    ((2*go.nX-3)*(1+go.alpha)+go.Ratio^2*(2*go.nY-3)*(1+go.beta));
go.aY = sqrt(kr2*go.sm2);
go.aX = go.Ratio*go.aY;

if minn > 2  % Fourth Spectral Moment
    CX4 = (3/2)*(go.nX-2)/((2*go.nX-3)^2*(2*go.nX-5));
    CY4 = (3/2)*(go.nY-2)/((2*go.nY-3)^2*(2*go.nY-5));
    CXY4 = (3/2)/((2*go.nX-3)*(2*go.nY-3));
    go.sm4 = CX*CX4*go.aX^4 + CY*CY4*go.aY^4 + ...
        CXY*(CX4*go.aX^4 + CY4*go.aY^4 + CXY4*go.aX^2*go.aY^2);
end

if nargin == 2
    psd = CX*bpqpsdx(go.nX,go.aX,f) + CY*bpqpsdx(go.nY,go.aY,f) + ...
        CXY*bpqpsdw(go.nX,go.nY,go.aX,go.aY,f);
end

end

Contact us