对于一个给定的随机粗糙面,对光波来说可能呈现很粗糙,而对微波来说却可能呈现的很光滑,这主要是因为随机表面的粗糙度是以波长为度量单位的统计参数来表征的。而描述粗糙面的统计量除功率谱密度外,还有高度起伏的概率密度函数和均方高,相关函数和相关长度,结构函数,特征函数,均方斜度与曲率半径等。下面将逐一介绍它们统计描述,并说明它们和功率谱密度间的相互关系。
图1 一维粗糙表面的模拟图(其中x轴和y轴都是以入射波长
1 粗糙面中各个基本参量
图1为一个一维粗糙表面,其高度起伏
粗糙面上任意两点间的关联程度,称为相关函数
其中
对表面高度相关函数作傅立叶变换,即得到其功率谱密度函数,表示为
对高斯相关函数和指数相关函数分别作如上傅立叶变换,得到高斯谱密度函数和指数谱密度函数,分别如下
均方根斜率
结构函数定义为表面上两点高度差的均值,即
上面费了好长时间,公式神马的实在是太麻烦了,真心想了解的同学可以下载下面的附件,我在这里就直接贴出代码了
http://files.cnblogs.com/xd-jinjian/%E7%AC%AC%E4%BA%8C%E7%AB%A0%EF%BC%88%E5%BB%BA%E6%A8%A1%EF%BC%89.pdf
第一种方法 不使用傅里叶变换,直接叠加法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102 |
// OneDimensionRoughSurface.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include <iostream> #include <complex> #include <stdlib.h> #include <fstream> #include <time.h> using
namespace std; const
complex< double > unit_i(0.0,1.0); const
complex< double > unit_r(1.0,0.0); const
double wavespeed=3.0* pow (10.0,8); //波速 const
double frequence=1.0* pow (10.0,9); //频率 const
double lamda=wavespeed/frequence; //波长 const
double delta=0.2*lamda; //均方根高度 const
double lcor=lamda; //相关长度 const
int M=1024; //总采样点数 const
double m=8; //一个波长内的采样点数 const
double dx=lamda/m; //x方向的间隔 const
double L=(M-1)*dx; //const double L=102.4*lamda; const
double pi=3.1415926; void
random( double
start, double
end, double
seed, double
size, double
*myrand); double
gaussrand(); int
_tmain( int
argc, _TCHAR* argv[]) { double
*x= new
double [M]; for ( int
i=0;i<M;++i) x[i]=(-0.5*M+i)*dx; double
*kj= new
double [M]; complex< double > *Fkj= new
complex< double > [M]; double
*S= new
double [M]; for ( int
i=0;i<M;++i) { kj[i]=(-0.5*M+i)*2*pi/L; S[i]=delta*delta*lcor/(2* sqrt (pi))* exp (-kj[i]*kj[i]*lcor*lcor/4); } double
*myrand= new
double [M]; random(0.0,1.0, time (0),M,myrand); Fkj[0]= sqrt (2*pi*L*S[0])*myrand[0]; Fkj[M/2]= sqrt (2*pi*L*S[M/2])*myrand[1]; for ( int
i=1;i<M/2;++i) { Fkj[i]= sqrt (2*pi*L*S[i])*(myrand[2*i]+unit_i*myrand[2*i+1])/ sqrt (2.0); Fkj[M-i]=conj(Fkj[i]); } complex< double > *res= new
complex< double > [M]; for ( int
i=0;i<M;++i) { double
xn=(-0.5*M+i)*dx; //double xn=(i)*dx; for ( int
j=0;j<M;++j) { res[i]+=Fkj[j]* exp (unit_i*kj[j]*xn); } res[i]/=L; } ofstream sout( "diejia.dat" ,ios::out); for ( int
i=0;i<M;++i) sout<<res[i].real()/lamda<< " " <<res[i].imag()/lamda<<endl; sout.close(); delete
[] res; delete
[] x; delete
[] S; delete
[] myrand; delete
[] kj; delete
[] Fkj; return
0; } void
random( double
start, double
end, double
seed, double
size, double
*myrand) //在用二维数组作为形参时,必须指明列数 { double
s=65536.0; double
w=2053.0; double
v=13849.0; double
T=0.0; int
m; for ( int
i=0;i<size;i++) { T=0.0; for ( int
k=0;k<12;k++) { seed=w*seed+v; m=seed/s; seed=seed-m*s; T=T+seed/s; } myrand[i]=start+end*(T-6.0); } } |
第二种,使用傅里叶变换
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176 |
// OneDimensionRoughSurface.cpp : 定义控制台应用程序的入口点。 // #include "stdafx.h" #include <iostream> #include <complex> #include <stdlib.h> #include <fstream> #include <time.h> using
namespace std; const
complex< double > unit_i(0.0,1.0); const
complex< double > unit_r(1.0,0.0); const
double wavespeed=3.0* pow (10.0,8); //波速 const
double frequence=1.0* pow (10.0,9); //频率 const
double lamda=wavespeed/frequence; //波长 const
double delta=0.2*lamda; //均方根高度 const
double lcor=lamda; //相关长度 const
int M=1024; //总采样点数 const
double m=8; //一个波长内的采样点数 const
double dx=lamda/m; //x方向的间隔 const
double L=(M-1)*dx; const
double pi=3.1415926; void
random( double
start, double
end, double
seed, double
size, double
*myrand); double
gaussrand(); inline
void swap ( float
&a, float
&b); void
bitrp ( float
xreal [], float
ximag [], int
n); void
IFFT ( float
xreal [], float
ximag [], int
n); int
_tmain( int
argc, _TCHAR* argv[]) { double
*x= new
double [M]; for ( int
i=0;i<M;++i) x[i]=-L/2+i*dx; double
*kj= new
double [M]; complex< double > *Fkj= new
complex< double > [M]; double
*S= new
double [M]; for ( int
i=0;i<=M/2;++i) { kj[i]=i*2*pi/L; S[i]=delta*delta*lcor/(2* sqrt (pi))* exp (-kj[i]*kj[i]*lcor*lcor/4); } double
*myrand= new
double [M]; random(0.0,1.0, time (0),M,myrand); Fkj[0]= sqrt (2*pi*L*S[0])*myrand[0]; Fkj[M/2]= sqrt (2*pi*L*S[M/2])*myrand[1]; for ( int
i=1;i<M/2;++i) { Fkj[i]= sqrt (2*pi*L*S[i])*(myrand[2*i]+unit_i*myrand[2*i+1])/ sqrt (2.0); Fkj[M-i]=conj(Fkj[i]); } float
xreal [M] = {}, ximag [M] = {}; for ( int
i=0;i<M;++i) { xreal[i]=Fkj[i].real(); ximag[i]=Fkj[i].imag(); } IFFT (xreal, ximag, M); ofstream sout( "fft.dat" ,ios::out); for ( int
i=0;i<M;++i) { sout<<xreal[i]*M/L/lamda<< " " <<ximag[i]*M/L/lamda<<endl; } sout.close(); delete
[] x; delete
[] S; delete
[] myrand; delete
[] kj; delete
[] Fkj; return
0; } void
random( double
start, double
end, double
seed, double
size, double
*myrand) //在用二维数组作为形参时,必须指明列数 { double
s=65536.0; double
w=2053.0; double
v=13849.0; double
T=0.0; int
m; for ( int
i=0;i<size;i++) { T=0.0; for ( int
k=0;k<12;k++) { seed=w*seed+v; m=seed/s; seed=seed-m*s; T=T+seed/s; } myrand[i]=start+end*(T-6.0); } } inline
void swap ( float
&a, float
&b) { float
t; t = a; a = b; b = t; } void
bitrp ( float
xreal [], float
ximag [], int
n) { // 位反转置换 Bit-reversal Permutation int
i, j, a, b, p; for
(i = 1, p = 0; i < n; i *= 2) { p ++; } for
(i = 0; i < n; i ++) { a = i; b = 0; for
(j = 0; j < p; j ++) { b = (b << 1) + (a & 1); // b = b * 2 + a % 2; a >>= 1; // a = a / 2; } if
( b > i) { swap (xreal [i], xreal [b]); swap (ximag [i], ximag [b]); } } } void
IFFT ( float
xreal [], float
ximag [], int
n) { // 快速傅立叶逆变换 float
wreal [M / 2], wimag [M / 2], treal, timag, ureal, uimag, arg; int
m, k, j, t, index1, index2; bitrp (xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 arg = 2 * pi / n; treal = cos
(arg); timag = sin
(arg); wreal [0] = 1.0; wimag [0] = 0.0; for
(j = 1; j < n / 2; j ++) { wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; } for
(m = 2; m <= n; m *= 2) { for
(k = 0; k < n; k += m) { for
(j = 0; j < m / 2; j ++) { index1 = k + j; index2 = index1 + m / 2; t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2]; timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; ureal = xreal [index1]; uimag = ximag [index1]; xreal [index1] = ureal + treal; ximag [index1] = uimag + timag; xreal [index2] = ureal - treal; ximag [index2] = uimag - timag; } } } for
(j=0; j < n; j ++) { xreal [j] /= n; ximag [j] /= n; } } |
结果如上图所示
原文:http://www.cnblogs.com/xd-jinjian/p/3674149.html