Skip to content
Snippets Groups Projects
Commit 721a20e0 authored by Jacobs Matthieu's avatar Jacobs Matthieu
Browse files

Matlab Scripts for artificial noise

parent bbc84030
No related branches found
No related tags found
No related merge requests found
% 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(p,r);
% 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;
[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');
% 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;
[P,R,T] = meshgrid(p,r,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)*nr + (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);
% To extract the noise, store in same format as BFSTREN (1 noise value for
% every grid point) => then shift to cell centers
% 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
% 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)
% 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
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
% 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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment