diff --git a/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.asv b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.asv new file mode 100644 index 0000000000000000000000000000000000000000..704454778a0c369e9340b31e6f685f1fe82ae17b --- /dev/null +++ b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid.asv @@ -0,0 +1,119 @@ +% Generate noise with spatial correlation in poloidal, radial and toroidal +% directions +% Sampling on spherical grid with spatial correlation +close all +clear all +clc + +Rmax = 500; +rmin = 100; +rstart = 10; + +L = 2*pi*Rmax; +% Parameters of the noise distribution +mu = 0; % mean of the noise (spatially constant) +sigma = 1; % standard deviation for the noise (spatially constant) +Lr = rmin/5; % correlation length in radial direction +Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction +Lt = 2*pi*Rmax/2; % Correlation "length" in toroidal direction + + +nr = 10; % number of grid points in x-direction +np = 5; % number of grid points in y-direction +nt = 10; + +r = (rstart/rmin:0.9/(nr-1):1)*rmin; +p = (0:1/(np-1):1)*2*pi; +[P,R] = meshgrid(r,p); + +% Vector of mean at each grid point +Mu = mu*ones(nr*np,1); + +% Construct covariance matrices +Sigma_corr = zeros(nr*np); +Sigma_no_corr = zeros(nr*np); + +for j = 1:np + for i = 1:nr + + k = i + (j-1)*nr; + + % Main diagonal: variance + Sigma_no_corr(k,k) = sigma^2; + Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + for n = 1:np + for m = 1:nr + l = m + (n-1)*nr; + if l ~= k + Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n)-R(i,j))^2)/Lr)*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp); + end + end + end + end +end + +b = issymmetric(Sigma_corr) +d = eig(Sigma_corr) +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nr,np); +noise_corr = reshape(noise_corr,nr,np); + +% Plot +figure;surf(R,P,noise_no_corr);title('Independent'); +figure;surf(R,P,noise_corr); title('Spatial correlation'); + +% Also check toroidal + +t = (0:1/(nt-1):1)*2*pi; +p = (0:1/(np-1):1)*2*pi; +[T,P] = meshgrid(t,p); + +% Vector of mean at each grid point +Mu = mu*ones(nt*np,1); + +% Construct covariance matrices +Sigma_corr = zeros(nt*np); +Sigma_no_corr = zeros(nt*np); + +for j = 1:np + for i = 1:nt + + k = i + (j-1)*nt; + + % Main diagonal: variance + Sigma_no_corr(k,k) = sigma^2; + Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + for n = 1:np + for m = 1:nt + l = m + (n-1)*np; + if l ~= k + Sigma_corr(k,l) = sigma^2*exp(-(abs(T(m,n)-T(i,j))*(Rmax+R(m,n)*cos(P(m,n))+Rmax+R(i,j)*cos(P(i,j)))/2)/Lt) ... + *exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp); + end + end + end + end +end + +b = issymmetric(Sigma_corr) +d = eig(Sigma_corr) +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nt,np); +noise_corr = reshape(noise_corr,nt,np); + +% Plot +figure;surf(T,P,noise_no_corr);title('Independent'); +figure;surf(T,P,noise_corr); title('Spatial correlation'); + diff --git a/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.asv b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.asv new file mode 100644 index 0000000000000000000000000000000000000000..8dfb478f69a64f5e9fe992fb2d421f31ee106dd5 --- /dev/null +++ b/Matlab Noise Analysis/CorrelatedNoiseSphericalGrid3D.asv @@ -0,0 +1,77 @@ +% Generate noise with spatial correlation in poloidal, radial and toroidal +% directions +% Sampling on spherical grid with spatial correlation +% 3D case +close all +clear all +clc + +Rmax = 500; +rmin = 100; +rstart = 10; + +L = 2*pi*Rmax; +% Parameters of the noise distribution +mu = 0; % mean of the noise (spatially constant) +sigma = 1; % standard deviation for the noise (spatially constant) +Lr = rmin/10; % correlation length in radial direction +Lp = 2*pi*rmin/10; % Correlation "length" in poloidal direction +Lt = 2*pi*Rmax/10; % Correlation "length" in toroidal direction + + +nr = 5; +np = 10; +nt = 2; + +r = (rstart/rmin:0.9/(nr-1):1)*rmin; +p = (0:1/(np-1):1)*2*pi; +t = (0:1/(nt-1):1)*2*pi; +[R,P,T] = meshgrid(r,p,t); + +% Vector of mean at each grid point +Mu = mu*ones(nr*np*nt,1); + +% Construct covariance matrices +Sigma_corr = zeros(nr*np*nt); +Sigma_no_corr = zeros(nr*np*nt); + +for h = 1:nt + for j = 1:np + for i = 1:nr + + k = i + (j-1)*nr + (h-1)*nr*np; + + % Main diagonal: variance + Sigma_no_corr(k,k) = sigma^2; + Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + for p = 1:nt + for n = 1:np + for m = 1:nr + l = m + (n-1)*np + (p-1)*nr*np; + if l ~= k + Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n,p)-R(i,j,h))^2)/Lr)... + *exp(-(abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2)/Lp) ... + *exp(-(abs(T(m,n,p)-T(i,j,h))*(Rmax+R(m,n,p)*cos(P(m,n,p))+Rmax+R(i,j,h)*cos(P(i,j,h)))/2)/Lt); + end + end + end + end + end + end +end + +b = issymmetric(Sigma_corr) +d = eig(Sigma_corr) +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nr,np,nt); +noise_corr = reshape(noise_corr,nr,np,nt); + +% Plot +%figure;surf(R,P,T);title('Independent'); +%figure;surf(R,P,T); title('Spatial correlation'); \ No newline at end of file diff --git a/Matlab Noise Analysis/Sparse2DWouter.m b/Matlab Noise Analysis/Sparse2DWouter.m new file mode 100644 index 0000000000000000000000000000000000000000..436cbb7e512d4a807f301ae9aef0986c77573e27 --- /dev/null +++ b/Matlab Noise Analysis/Sparse2DWouter.m @@ -0,0 +1,65 @@ +% Sampling on cartesian grid with spatial correlation +close all +clear all +clc + + +% Parameters of the noise distribution +mu = 0; % mean of the noise (spatially constant) +sigma = 1; % standard deviation for the noise (spatially constant) +L = 0.25; % correlation length (spatially constant) + +Ln =0.65; +% Grid to sample on +Lx = 1; % length of domain in x-direction +Ly = 1; % length of domain in y-direction + +nx = 10; % number of grid points in x-direction +ny = 10; % number of grid points in y-direction + +x = (0:1/(nx-1):1)*Lx; +y = (0:1/(ny-1):1)*Ly; +[X,Y] = meshgrid(x,y); + + +% Vector of mean at each grid point +Mu = mu*ones(nx*ny,1); + +% Construct covariance matrices +Sigma_corr = zeros(nx*ny); +Sigma_no_corr = zeros(nx*ny); + +for j = 1:ny + for i = 1:nx + + k = i + (j-1)*ny; + + % Main diagonal: variance + Sigma_no_corr(k,k) = sigma^2; + Sigma_corr(k,k) = sigma^2; + % Off-diagonal: spatial correlation + for n = 1:ny + for m = 1:nx + l = m + (n-1)*ny; + if l ~= k + if ((X(m,n)-X(i,j))^2 + (Y(m,n)-Y(i,j))^2)^0.5 < Ln + Sigma_corr(k,l) = sigma^2*exp(-((X(m,n)-X(i,j))^2 + (Y(m,n)-Y(i,j))^2)^0.5/L); + end + end + end + end + end +end + + +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nx,ny); +noise_corr = reshape(noise_corr,nx,ny); + +% Plot +figure;surf(X,Y,noise_no_corr);title('Independent'); +figure;surf(X,Y,noise_corr); title('Spatial correlation'); \ No newline at end of file diff --git a/Matlab Noise Analysis/SparseCorr2D.m b/Matlab Noise Analysis/SparseCorr2D.m new file mode 100644 index 0000000000000000000000000000000000000000..e6a98255517cd3a444a2553d994bda867cd0e169 --- /dev/null +++ b/Matlab Noise Analysis/SparseCorr2D.m @@ -0,0 +1,146 @@ +% Generate noise with spatial correlation in poloidal, radial and toroidal +% directions +% Sampling on spherical grid with spatial correlation +close all +clear all +clc + +Rmax = 500; +rmin = 100; +rstart = 10; + +L = 2*pi*Rmax; +% Parameters of the noise distribution +mu = 0; % mean of the noise (spatially constant) +sigma = 1; % standard deviation for the noise (spatially constant) +Lr = rmin/10; % correlation length in radial direction +Lp = 2*pi*rmin/20; % Correlation "length" in poloidal direction +Lt = 2*pi*Rmax; % Correlation "length" in toroidal direction +%Lp = Lr; +Lcorr = Lr/2; + +nr = 100; % number of grid points in x-direction +np = 100; % number of grid points in y-direction +nt = 10; + +r = (rstart/rmin:0.9/(nr-1):1)*rmin; +p = (0:1/(np-1):1)*2*pi; +[P,R] = meshgrid(p,r); + +% Vector of mean at each grid point +Mu = mu*ones(nr*np,1); + +is = (1:np*nr); +js = (1:np*nr); +vs = sigma^2*ones(np*nr,1); + +% Construct covariance matrices +Sigma_corr = sparse(is,js,vs,nr*np,nr*np); +%Sigma_corr = zeros(nr*np); +Sigma_no_corr = sparse(is,js,vs,nr*np,nr*np); + +for j = 1:np + for i = 1:nr + + k = i + (j-1)*nr; + + % Main diagonal: variance + % Sigma_no_corr(k,k) = sigma^2; + % Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + nnb = 4; + % for n = j-nnb:j+nnb + % for m = i-nnb:i+nnb + % if (m ==i || n==j) && not(m==i && n==j) %&& not(m==0 || m==nr+1 ||n==0 ||n == np+1) + % mc = m+nr*(m<=0) - nr*(m>=nr+1); + % nc = n+np*(n<=0) - np*(n>=np+1); + % l = mc + (nc-1)*nr; + % Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(mc,nc)-R(i,j))^2)/Lr)*exp(-(abs(P(mc,nc)-P(i,j))*(R(mc,nc)+R(i,j))/2)/Lp); + % end + % end + %end + % Alternatief + for n=1:np + for m=1:nr + l = m+(n-1)*nr; + if l ~= k + distance = ((R(m,n)-R(i,j))^2+ (abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2))^0.5; + if distance < Lcorr + Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n)-R(i,j))^2)/Lr)*exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp); + %Sigma_corr(k,l) = sigma^2*exp(-distance/Lcorr); + end + end + end + end + + end +end + +b = issymmetric(Sigma_corr) +%d = eig(Sigma_corr) + +[U,S,V] = svd(Sigma_corr); +z = 0+randn(nr*np,1)*sigma; +noise_corr = U *sqrt(S)*z; +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +%noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nr,np); +noise_corr = reshape(noise_corr,nr,np); + +% Plot +figure;surf(R,P,noise_no_corr);title('Independent'); +figure;surf(R,P,noise_corr); title('Spatial correlation'); + +% Also check toroidal + +t = (0:1/(nt-1):1)*2*pi; +p = (0:1/(np-1):1)*2*pi; +[P,T] = meshgrid(p,t); + +% Vector of mean at each grid point +Mu = mu*ones(nt*np,1); + +% Construct covariance matrices +Sigma_corr = zeros(nt*np); +Sigma_no_corr = zeros(nt*np); + +for j = 1:np + for i = 1:nt + + k = i + (j-1)*nt; + + % Main diagonal: variance + Sigma_no_corr(k,k) = sigma^2; + Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + for n = 1:np + for m = 1:nt + l = m + (n-1)*nt; + if l ~= k + Sigma_corr(k,l) = sigma^2*exp(-(abs(T(m,n)-T(i,j))*(Rmax+R(m,n)*cos(P(m,n))+Rmax+R(i,j)*cos(P(i,j)))/2)/Lt) ... + *exp(-(abs(P(m,n)-P(i,j))*(R(m,n)+R(i,j))/2)/Lp); + end + end + end + end +end + +b = issymmetric(Sigma_corr) +d = eig(Sigma_corr) +% Generate the noise from multivariate normal distributions +noise_no_corr = mvnrnd(Mu,Sigma_no_corr); +noise_corr = mvnrnd(Mu,Sigma_corr); + +% Reshape onto mesh +noise_no_corr = reshape(noise_no_corr,nt,np); +noise_corr = reshape(noise_corr,nt,np); + +% Plot +figure;surf(T,P,noise_no_corr);title('Independent'); +figure;surf(T,P,noise_corr); title('Spatial correlation'); + diff --git a/Matlab Noise Analysis/SparseCorr3D.m b/Matlab Noise Analysis/SparseCorr3D.m new file mode 100644 index 0000000000000000000000000000000000000000..238a1e732288cb7d5243594b4898e132d894fa9f --- /dev/null +++ b/Matlab Noise Analysis/SparseCorr3D.m @@ -0,0 +1,105 @@ +% Generate noise with spatial correlation in poloidal, radial and toroidal +% directions +% Sampling on spherical grid with spatial correlation +close all +clear all +clc + +Rmax = 500; +rmin = 100; +rstart = 10; + +L = 2*pi*Rmax; +% Parameters of the noise distribution +mu = 0; % mean of the noise (spatially constant) +sigma = 1; % standard deviation for the noise (spatially constant) +Lr = rmin/10; % correlation length in radial direction +Lp = 2*pi*rmin/20; % Correlation "length" in poloidal direction +Lt = 2*pi*Rmax; % Correlation "length" in toroidal direction +%Lp = Lr; +Lcorr = Lr/2; + +nr = 10; % number of grid points in x-direction +np = 10; % number of grid points in y-direction +nt = 10; +tol = 10^-15; +r = (rstart/rmin:0.9/(nr-1):1)*rmin; +p = (0:1/(np-1):1)*2*pi; +t = (0:1/(nt-1):1)*2*pi; +[P,R,T] = meshgrid(p,r,t); +% Vector of mean at each grid point +Mu = mu*ones(nr*np*nt,1); + +is = (1:np*nr*nt); +js = (1:np*nr*nt); +vs = sigma^2*ones(np*nr*nt,1); + +% Construct covariance matrices +Sigma_corr = sparse(is,js,vs,nr*np*nt,nr*np*nt); +%Sigma_corr = zeros(nr*np); +Sigma_no_corr = sparse(is,js,vs,nr*np*nt,nr*np*nt); + +for h = 1:nt + for j = 1:np + for i = 1:nr + + k = i + (j-1)*nr + (h-1)*nr*np; + + % Main diagonal: variance + % Sigma_no_corr(k,k) = sigma^2; + % Sigma_corr(k,k) = sigma^2; + + % Off-diagonal: spatial correlation + nnb = 4; + % for n = j-nnb:j+nnb + % for m = i-nnb:i+nnb + % if (m ==i || n==j) && not(m==i && n==j) %&& not(m==0 || m==nr+1 ||n==0 ||n == np+1) + % mc = m+nr*(m<=0) - nr*(m>=nr+1); + % nc = n+np*(n<=0) - np*(n>=np+1); + % l = mc + (nc-1)*nr; + % Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(mc,nc)-R(i,j))^2)/Lr)*exp(-(abs(P(mc,nc)-P(i,j))*(R(mc,nc)+R(i,j))/2)/Lp); + % end + % end + %end + % Alternatief + for p = 1:nt + for n=1:np + for m=1:nr + l = m+(n-1)*nr+(p-1)*nr*np; + if l ~= k + distance = ((R(m,n,p)-R(i,j,h))^2+ (abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2))^0.5; + if distance < Lcorr + Sigma_corr(k,l) = sigma^2*exp(-sqrt((R(m,n,p)-R(i,j,h))^2)/Lr)... + *exp(-(abs(P(m,n,p)-P(i,j,h))*(R(m,n,p)+R(i,j,h))/2)/Lp)... + *exp(-(abs(T(m,n,p)-T(i,j,h))*(Rmax+R(m,n,p)*cos(P(m,n,p))+Rmax+R(i,j,h)*cos(P(i,j,h)))/2)/Lt); + b=1; % test bp + %Sigma_corr(k,l) = sigma^2*exp(-distance/Lcorr); + end + end + end + end + end + + end + end +end + + b = issymmetric(Sigma_corr) + %d = eig(Sigma_corr) + + [U,S,V] = svds(Sigma_corr,nr*np*nt); + z = 0+randn(nr*np*nt,1)*sigma; + noise_corr = U *sqrt(S)*z; + % Generate the noise from multivariate normal distributions + noise_no_corr = mvnrnd(Mu,Sigma_no_corr); + %noise_corr = mvnrnd(Mu,Sigma_corr); + + % Reshape onto mesh + noise_no_corr = reshape(noise_no_corr,nr,np,nt); + noise_corr = reshape(noise_corr,nr,np,nt); + + % Plot + %figure;surf(R,P,noise_no_corr);title('Independent'); + %figure;surf(R,P,noise_corr); title('Spatial correlation'); + + \ No newline at end of file