Posted Updated Codes3 minutes read (About 406 words)

多径衰落仿真

做实验时写出来的代码，感觉会用的上，记录一下。

task1 - 计算多径衰落的时域分布。

task2 - 计算多径衰落的空间分布。

task3 - 用BPSK调制，对比有无衰落的误码率。

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
clear;
close all;
clc;

L = 1:1000;
N = 16;
v = 40 / 3.6;
f_c = 800e6;
C = physconst('lightspeed');
R = 9600;

ramda = C / f_c;
T = 1 / R;
FDT = v * f_c / C * T;

% --- Task 1

theta = zeros([N, 1]);
for n = 1 : N
    if(n == 1)
        theta(n) = rand() * 2 * pi / 16;
    else
        theta(n) = theta(n - 1) + rand() * 2 * pi / 16;
    end
end

phi = rand([N, 1]) * 2 * pi;

alpha = 2 * pi * FDT * cos(theta);

k = 1 : 1000;

Fk = exp(1j * ((alpha * k) + phi));
Fk = sum(Fk, 1) / sqrt(N);
Fk = 20 * log10(abs(Fk));

figure(1);
plot(Fk);
xlabel('k'), ylabel('F(k)/dB');

% --- Task 2

X = 1 : 100;
Y = 1 : 100;

F = zeros(100, 100);
for x = X
    for y = Y
        F(x, y) = abs(sum(exp(1j * ((y / 100 * cos(theta) - x / 100 * sin(theta)) * 2 * pi / ramda + phi))) / sqrt(N));
    end
end

figure(2);
mesh(F);
xlabel('x'), ylabel('y'), zlabel('F(x, y)');

% --- Task 3

snr_in_db = 0 : 20;
M = 1e6;

data = zeros([M, 1]);
for i = 1 : M
    temp = rand();
    if(temp > 0.5)
        data(i) = 1;
    else 
        data(i) = 0;
    end
end

signal = zeros([M, 1]);

for i = 1 : M
    signal(i) = 2 * data(i) - 1;
end

pf = zeros([length(snr_in_db), 1]);
for i = 1 : length(snr_in_db)
    numoferr = 0;
    
    k = 1 : M;
    
    Fk = exp(1j * ((alpha * k) + phi));
    Fk = sum(Fk, 1) / sqrt(N);
    
    rk = awgn(signal .* Fk', snr_in_db(i));
    
    dk = real(rk .* conj(Fk'));
    
    for k = 1 : M
        if (dk(k) > 0)
            a_cd = 1;
        else
            a_cd = 0;
        end
        
        if (a_cd ~= data(k))
            numoferr = numoferr + 1;
        end
    end
    
    pf(i) = numoferr / M;
end

pa = zeros([length(snr_in_db), 1]);
for i = 1 : length(snr_in_db)
    numoferr = 0;
    
    snr = power(10, snr_in_db(i) / 10);

    rk = signal + wgn(length(signal), 1, 10 * log10(1 / 2 / snr));
    
    dk = rk;
    
    for k = 1 : M
        if (dk(k) > 0)
            a_cd = 1;
        else
            a_cd = 0;
        end
        
        if (a_cd ~= data(k))
            numoferr = numoferr + 1;
        end
    end
    
    pa(i) = numoferr / M;
end

figure(3);
semilogy(snr_in_db, [pf, pa]);
grid on;
xlabel('SNR'), ylabel('BER');

运行结果

F(k)

F(x, y)

误码率对比（蓝：衰落 橙：无衰落）

Posted Updated Codes3 minutes read (About 444 words)

矩形平面阵阵因子计算

项目需求，需要算一下平面阵的阵因子。

本文计算阵因子使用的是比较简单的方法，利用直线阵的公式计算平面阵。

实际上还可以使用向量进行分析，更加的泛用。

均匀直线阵的阵因子公式如下：

$$ f(\psi) = \frac{sin(\frac{N}{2}\psi)}{sin(\frac{1}{2}\psi)} $$

$$ \psi = kdcos\delta + \xi $$

$d$是阵元间距，$\delta$是观察方向和阵轴的夹角，$\xi$是阵元的相移。

对于平面阵，我们可以把它分解为X方向和Y方向的两个均匀直线阵。由乘法原理可知：

$$ f(\theta, \varphi) = f_{1}(\theta, \varphi) \times f_{2x} (\theta, \varphi) \times f_{2y}(\theta, \varphi) $$

计算时只需要换算一下阵轴夹角即可：

$$ cos\delta_{x} = sin\theta sin\varphi $$

$$ cos\delta_{y} = sin\theta cos\varphi $$

代码如下

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
clear;
close all;
clc;

% 指定频率
freq = 1.5e6;

c = physconst('lightspeed');
lambda = c / freq;
k = 2 * pi / lambda;

% 指定平面阵参数
dx = 1 / 4 * lambda;
dy = 1 / 4 * lambda;
nx = 3;
ny = 3;
pdx = 0;
pdy = 0;

theta = meshgrid(eps : pi / 180 : pi);
phi = meshgrid(eps : 2 * pi / 180 : 2 * pi)';

f2x = sin(nx / 2 * (k * dx * (sin(theta) .* cos(phi)) + pdx)) ./ sin(1 / 2 * (k * dx * (sin(theta) .* cos(phi)) + pdx));
f2y = sin(ny / 2 * (k * dy * (sin(theta) .* sin(phi)) + pdy)) ./ sin(1 / 2 * (k * dy * (sin(theta) .* sin(phi)) + pdy));

f = f2x .* f2y;

f = abs(f);

[x, y, z] = sph2cart(phi, pi / 2 - theta, f);

mesh(x, y, z);

title('平面阵方向图');
xlabel('x'), ylabel('y'), zlabel('z');
axis equal;

由于Matlab不会去计算极限，所以计算时会出现NaN，无法计算对数坐标。

3X3阵

端射阵

网格上的小缺口就是NaN的位置，公式计算几乎无法消除这个问题。

之后我会尝试使用向量计算，这样就可以规避掉上面的问题。

Posted Updated Dailya minute read (About 114 words)

天线方向图

之前做实验，老师要求用Origin画天线的方向图，结果我实在是整不到破解版，就自己用Matalb整了一个。

代码就不放了，毕竟确实也比较简单。

效果还是不错的。

天线方向图

这里半功率主瓣也标出来了，可以用鼠标点一下上下的点来计算主瓣宽度。

Posted Updated Codes8 minutes read (About 1177 words)

Matlab史密斯圆图与LNA设计

这几天需要用Matlab设计一个LNA，结果发现Matlab的史密斯圆图不太好用…

于是我决定自己整一个史密斯圆图。

电阻圆族和电抗圆族公式如下：

$${(\Gamma_{r}-\frac{r}{1+r})}^2 + {\Gamma_{i}}^2 = {(\frac{1}{1+r})}^2$$

$${(\Gamma_{r}-1)}^2 + {(\Gamma_{i}-\frac{1}{x})}^2 = {(\frac{1}{x})}^2$$

那么，可以写出代码如下：

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
function [] = smith()
    for x = -4 : 0.2 : 4
        viscircles([1, 1 / x], 1 / abs(x), "LineWidth", 0.1, "Color", "blue");
        xlim([-1, 1]);
        ylim([-1, 1]);
    end
    line([-1, 1], [0, 0], "LineWidth", 0.1, "Color", "blue");
    
    for r = 0 : 0.2 : 5
        viscircles([r / (1 + r), 0], 1 / (1 + r), "LineWidth", 0.1, "Color", "blue");
        hold on; 
        axis equal;
        xlim([-1, 1]);
        ylim([-1, 1]);
    end     

    anglex={-1, -1, 1, 1};
    angley={1, -1, -1, 1};
    for i = 1 : 4 
        for theta = pi / 2 * i : pi / 1000 : pi / 2 * ( i + 1 )
            x = cos(theta); 
            y = sin(theta);
            x1 = cos(theta + pi / 1000); 
            y1 = sin(theta + pi / 1000);
            fill([x, x1, anglex{i}], [y, y1, angley{i}], "b");
        end
    end
end

效果是这样的

史密斯圆图

LNA设计基于上文的史密斯圆图来表示。

具体公式比较复杂，可以去看《射频电路设计――理论与应用（第二版）》这本书📕，公式位于第九章。

代码如下：

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
clear;
close all;
clc;

% S11 A11 S21 A21 S12 A12 S22 A22 Topt_m Topt_a Fmin_db Rn/50 顺序 
Parameter_demo = [0.3 30 2.5 -80 0.2 -60 0.2 -15 0.5 45 1.5 0.08];
Parameter_2G = [];

% 选择S参数
Sel = Parameter_demo;

% 指定 增益
G_db = 8;
G = power(10, G_db / 10);

% 指定 VSWRimn
VSWRimn = 1.5;

Timn_abs = (VSWRimn - 1) / (VSWRimn + 1);

% 指定 噪声系数
F_db = 1.6; 
F = power(10, F_db / 10);

% Rn/50
Rn = Sel(12);

% 最小噪声系数
Fmin_db = Sel(11); 
Fmin = power(10, Fmin_db / 10);

% opt反射系数
Topt_a = ang2rad(Sel(10));
Topt = Sel(9) * exp (1i * Topt_a);

% 指定S参数
a11 = ang2rad(Sel(2));
s11 = Sel(1) * exp (1i * a11);

a21 = ang2rad(Sel(4));
s21 = Sel(3) * exp (1i * a21);

a12 = ang2rad(Sel(6));
s12 = Sel(5) * exp (1i * a12);

a22 = ang2rad(Sel(8));
s22 = Sel(7) * exp (1i * a22);

% 史密斯圆图
subplot(121);
smith();

% 稳定性计算
delta = s11 * s22 -s12 * s21;

k = (1 - power(abs(s11), 2) - power(abs(s22), 2) + power(abs(delta), 2)) / (2 * abs(s12) * abs(s21));

% 最大功率增益
Gmax = abs(s21) / abs(s12) * (k  - sqrt(power(k, 2) - 1));
Gmax_db = 10 * log10(Gmax);

% 比例系数
g0 = G / power(abs(s21), 2);

% 功率增益圆
dg0 = g0 * conj(s22 - delta * conj(s11)) / (1 + g0 * (power(abs(s22), 2) - power(abs(delta), 2)));
rg0 = sqrt(1 - 2 * k * g0 * abs(s12 * s21) + power(g0, 2) * power(abs(s12 * s21), 2)) / abs(1 + g0 * (power(abs(s22), 2) - power(abs(delta), 2)));

% 转换
dgs = ((1 - s22 * dg0) * conj((s11 - delta * dg0)) - power(rg0, 2) * conj(delta) * s22) / (power(abs(1 - s22 * dg0), 2) - power(rg0, 2) * power(abs(s22), 2));
rgs = rg0 * abs(s12 * s21) / abs(power(abs(1 - s22 * dg0), 2) - power(rg0, 2) * power(abs(s22), 2));

viscircles([real(dgs), imag(dgs)], rgs, "LineWidth", 2, "Color", "red");

% 噪声功率圆
Qk = power(abs(1 + Topt), 2) * (F - Fmin) / (4 * Rn);
df = Topt / (1 + Qk);
rf = sqrt((1 - power(abs(Topt), 2)) * Qk + power(Qk, 2)) / (1 + Qk);

viscircles([real(df), imag(df)], rf, "LineWidth", 2, "Color", "green");

% 噪声opt
scatter(real(Topt), imag(Topt), "red");

% 扫描取优
min_F = inf;
VSWRomn_target = inf;
for degree = 1:360
    Ts = dgs + rgs * exp (1i * ang2rad(degree));
    Tl = (s11 - conj(Ts)) / (delta - s22 * conj(Ts));

    Tout = s22 + s12 * s21 * Ts / (1 - s11 * Ts);
    Tomn = (conj(Tout) - Tl) / (1 - Tl * Tout);

    VSWRomn = (1 + abs(Tomn)) / (1 - abs(Tomn));

    if(abs(Ts - df) < rf && VSWRomn < VSWRomn_target)
        scan_F = Fmin + 4 * Rn * power(abs(Ts - Topt), 2) / ((1 - power(abs(Ts), 2)) * power(abs(1 + Topt), 2));

        if(scan_F < min_F)
            min_F = scan_F;
            Ts_final = Ts;
            Tl_final = Tl;
        end
        scatter(real(Ts), imag(Ts), "yellow");
    end
end

scatter(real(Ts_final), imag(Ts_final), "red");

% 下面进行VSWR的优化

Tin = s11 + (s21 * s12 * Tl_final) / (1 - s22 * Tl_final);

dv_imn = (1 - power(Timn_abs, 2)) * conj(Tin) / (1 - power(Timn_abs * abs(Ts), 2));
rv_imn = (1 - power(abs(Tin), 2)) * Timn_abs / (1 - power(Timn_abs * abs(Ts), 2));

viscircles([real(dv_imn), imag(dv_imn)], rv_imn, "LineWidth", 2, "Color", "cyan");
VSWRomn_list = zeros(360, 1);

for degree = 1:360
    Ts = dv_imn + rv_imn * exp (1i * ang2rad(degree));

    Tout = s22 + s12 * s21 * Ts / (1 - s11 * Ts);
    Tomn = (conj(Tout) - Tl_final) / (1 - Tl_final * Tout);

    VSWRomn_list(degree) = (1 + abs(Tomn)) / (1 - abs(Tomn));
end

subplot(122);
grid on;
xlim([1, 360]);
ylim([1, inf]);
hold on;

% 画出不同角度的VSWR，然后选择一个角度
plot(1:360, VSWRomn_list, "Color", "blue");
line([1, 360], [VSWRimn, VSWRimn], "Color", "red", "LineStyle", "--");

prompt = "选择角度：";
degree = input(prompt);

Ts_final = dv_imn + rv_imn * exp (1i * ang2rad(degree));
Tout = s22 + s12 * s21 * Ts_final / (1 - s11 * Ts_final);

% Gt值计算
Gt = (1 - power(abs(Tl_final), 2)) * power(abs(s21), 2) * (1 - power(abs(Ts_final), 2)) / (power(abs(1 - Ts_final * Tin), 2) * power(abs(1 - s22 * Tl_final), 2));
Gt_db = 10 * log10(Gt);

% 噪声系数
F_final = Fmin + 4 * Rn * power(abs(Ts_final - Topt), 2) / ((1 - power(abs(Ts_final), 2)) * power(abs(1 + Topt), 2));
F_final_db = 10 * log10(F_final);

% 输出
fprintf("Ts: %f∠%f°\nTout: %f∠%f°\n", abs(Ts_final), getAngle(Ts_final), abs(Tout), getAngle(Tout));
fprintf("Gt(dB): %f\nF(dB): %f\nVSWRomn: %f\n", Gt_db, F_final_db, VSWRomn_list(degree));

该代码应配合下面的函数使用。

1
2
3
function [rad] = ang2rad(ang)
    rad = ang / 360 * 2 * pi;
end
1
2
3
function [ang] = getAngle(z)
    ang = angle(z) / 2 / pi * 360;
end

运行结果如下：

计算结果

左图中红圈就是对于$\Gamma_{S}$的等增益圆，绿圈则是等噪声圆。

绿圈中的红点是$\Gamma_{opt}$。

左图下方一大片黄点是满足参数要求（增益、噪声）的点，红色的则是噪声系数最小的点。

获得噪声系数最小的点后，进行驻波比优化。通过给定IMN的驻波比，画出图中蓝色的等驻波比圆。扫描该圆，计算驻波比，画出右图。

右图中实线是OMN的驻波比，需要找到一个比较小的点；虚线是IMN的驻波比。

将找到的OMN对应的角度输入Matlab，即可得到输出。

输出

Follow

Links

Categories

Recents

Archives

Tags

FPGA2
Game1
HLS1
HTTP-FLV1
LNA1
Matlab4
Photo1
QT2
RTMP1
Stable Diffusion1
TCP1
UART1
Unity1
Verilog2
分频器1
多径衰落1
数码管1
星穹铁道1
物联网1
矩形平面阵1
阵因子1
麦克斯韦方程组1
×