在Matlab中使用MVNRND創建2個高斯相關矩陣

Question

這是我在這裏詢問的第一個問題，請耐心等待。我對Matlab並不陌生，但之前從未使用過MVNRND函數，而且我的統計知識並不強大。我試圖做的總結如下：我試圖創建一個函數，生成2個相關的相位屏幕（NxN矩陣），將用於電磁高斯謝爾模型光束傳播模擬。光束需要用於X和Y偏振狀態的單獨的隨機相位屏幕。我到目前爲止的代碼如下。

function [phz_x,phz_y]=GSM_phase_screen_2(l_phi_x,l_phi_y,sigma_phi_x, ... 
sigma_phi_y,gamma,N,delta) 
%GSM_PHASE_SCREEN_2 
% This code generates two correlated 2-D NxN Gaussian Schell-Model (GSM) 
% phase screens (matrices) of unit variance circular complex Gaussian 
% random numbers for the X and Y polarization states provided l_phi_x, 
% l_phi_y, sigma_phi_x, sigma_phi_y, and gamma. It utilizes the MVNRND 
% Matlab function. 
% 
% l_phi_x:  [m] correlation length in X phase screen 
% l_phi_y:  [m] correlation length in Y phase screen 
% sigma_phi_x: phase variance in X phase screen 
% sigma_phi_y: phase variance in Y phase screen 
% gamma:   correlation coefficient 
% N:    number of samples per side of grid 
% delta:   [m] sample grid spacing 
% 
% phz_x:   [rad] 2-D phase screen for X polarization state 
% phz_y:   [rad] 2-D phase screen for Y polarization state 

% ORIGINAL AUTHOR: Santasri Basu 
% MODIFIED BY: Matthew Gridley 
%  Added input arguments needed to generate 2 correlated phase 
%  screens, updated PSD equations, and replaced RANDN with MVNRND 
% 

% Setup the Power Spectral Density (PSD) 
del_f = 1/(N*delta); % frequency grid spacing [1/m] 
fx = (-N/2 : N/2-1) * del_f; 

% Frequency grid 
[fx,fy] = meshgrid(fx); 

% GSM phase PSD 
PSD_phi_x = (sigma_phi_x^2) * pi * (l_phi_x^2) * gamma * ... 
    exp(-((pi * l_phi_x)^2) * (fx.^2 + fy.^2)); 
PSD_phi_y = (sigma_phi_y^2) * pi * (l_phi_y^2) * gamma * ... 
    exp(-((pi * l_phi_y)^2) * (fx.^2 + fy.^2)); 

% Random draws of Fourier series coefficients 
% (zero mean Gaussian random numbers) 
% 
% the 2 lines of code below need changed to generate the correlated random 
% draws using MVNRND and GAMMA 
cn_x = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi_x) * del_f; 
cn_y = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi_y) * del_f; 

% Synthesize the phase screens 
phz_x = real(ift2(cn_x,1)); 
phz_y = real(ift2(cn_y,1)); 

end 



function [g, x] = ift2(G, df) 
% [g, x] = ift2(G, df) 
% 2-D inverse Fourier transform that keeps the origin at the center of 
% the grid. 
% 
% G: Complex field in frequency space 
% df: Spacing in frequency space [m^-1] 
% g: Complex field in coordinate space 

% Core function written by Jason Schmidt 
% Modified: 17 Apr 2010 
% By: Daniel J. Wheeler 
% 
% x output added functionality by Michael Steinbock 6/8/2014 
% 

g = ifftshift(ifft2(ifftshift(G))) * (length(G) * df)^2; 

%% Calc x: 
if nargout == 2 
    N = size(G, 1); 

    x = (0 : N-1)/(N*delta_f); 
end 

在上面的代碼，代碼的兩行與開始的評論下面的「％隨機抽取傅立葉級數係數​​」是，我需要幫助。我以前用你看到的代碼做了兩個矩陣，但意識到它們不是高斯相關的。從我的學術顧問的建議，我應該使用MVNRND生成這些階段屏幕。在查看了MVNRND的幫助文件之後，我迷失在如何使用它的目的。我在這裏搜索，試圖找到類似的問題和答案，沒有運氣，我也搜索了谷歌。任何人都可以協助改變這兩行代碼來利用MVNRND。謝謝！

Answer 1

在這裏輸入代碼經過大量的研究，我想出瞭如何使用MVNRND來滿足我的目的。我的意圖是創建4個隨機NxN矩陣來替換下面的代碼片段中randn（N）的4個用法。我需要用MVNRND替換這些隨機矩陣的原因是它們是相關的。在MVNRND中，你必須提供一個協方差矩陣。這就是困擾我的東西。

cn_x = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi_x) * del_f; 
cn_y = (randn(N) + 1i*randn(N)) .* sqrt(PSD_phi_y) * del_f; 

要創建它，我有4個不同的隨機值，其中我必須計算組合對的方差（協方差）：rx_real，rx_imag，ry_real，和ry_imag。一旦我明白了這一點，我就可以創建協方差矩陣。

下一個問題是確定需要設置MVNRND中的'個案'值。我需要4個相關的NxN矩陣，所以我確定案例需要是4xN^2的矩陣。然後，我可以使用'reshape'命令將MVNRND輸出轉換爲需要的4個NxN矩陣。

請參閱下面的代碼。希望這可以幫助別人！

% Multivariate normal parameters 
mu = zeros([1,4]); % Zero mean Gaussian 
% Covariance matrix for 4 circular complex Gaussian random numbers: 
% rx_real, rx_imag, ry_real, ry_imag 
% 
% [<rx_real rx_real> <rx_real rx_imag> <rx_real ry_real> <rx_real ry_imag>; 
% <rx_imag rx_real> <rx_imag rx_imag> <rx_imag ry_real> <rx_imag ry_imag>; 
% <ry_real rx_real> <ry_real rx_imag> <ry_real ry_real> <ry_real ry_imag>; 
% <ry_imag rx_real> <ry_imag rx_imag> <ry_imag ry_real> <ry_imag ry_imag>] 
sigma = [1 0 gamma 0; 
     0 1 0 gamma; 
     gamma 0 1 0; 
     0 gamma 0 1]; 
cases = N^2; % matrix of random vectors 

r = mvnrnd(mu, sigma, cases); % gives a 512^2x4 double matrix 

rx_real = reshape(r(:,1),[N N]); 
rx_imag = reshape(r(:,2),[N N]); 
ry_real = reshape(r(:,3),[N N]); 
ry_imag = reshape(r(:,4),[N N]); 

% Correlated random draws of Fourier series coefficients 
cn_x = (rx_real + 1i*rx_imag) .* sqrt(PSD_phi_x) * del_f; 
cn_y = (ry_real + 1i*ry_imag) .* sqrt(PSD_phi_y) * del_f;