Posted Updated Codes2 minutes read (About 338 words)

UART传输协议实现

FPGA上的UART协议。

RX端的时钟是16×波特率；TX端则是等于波特率

本代码数据格式： 1 起始 8 数据 1 结束

RX端代码

clk - 时钟

rst - 重置 低有效

din - 接收信号线

srecv - 读取信号 上升沿触发

data - 接收数据

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
module uart_rx(
	input clk,
	input rst,
	input din,
	output reg srecv = 0,
	output reg [7:0] data = 0);
	
	reg [3:0] count = 0;
	reg [1:0] state = 0;
	reg [2:0] i = 0;
	
	always @ (posedge clk, negedge rst) begin
		if(~rst) begin
			count <= 0;
			i <= 0;
			srecv <= 0;
			data <= 0;
			state <= 2'h0;
		end
		else begin
			case(state)
				2'h0: begin
					if(din == 0)
						if(count != 15)
							count <= count + 1;
						else begin
							count <= 0;
							state <= 2'h1;
							srecv <= 0;
						end
					else
						count <= 0;
				end
				
				2'h1: begin
					if(count != 15) begin
						count <= count + 1;
						if(count == 10) begin
							data[i] <= din;
						end
					end
					else begin
						count <= 0;

						if(i == 7)
							state <= 2'h2;

						i <= i + 1;
					end
				end
				
				2'h2: begin
					if(count != 15)
						count <= count + 1;
					else begin
						count <= 0;
						state <= 2'h0;
						srecv <= 1;
					end
				end
				
			endcase
		end
	end
	
endmodule

TX端代码

clk - 时钟

rst - 重置 低有效

ssend - 发送信号 上升沿触发

data - 发送数据

send - 发送状态 高为正在发送

dout - 发送信号线

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
module uart_tx(
	input clk,
	input rst,
	input ssend,
	input [7:0] data,
	output reg send = 0,
	output reg dout = 1);
	
	reg [1:0] state = 2'h2;
	reg [7:0] sdata = 0;
	reg [3:0] i = 0;
	
	reg last_check = 0;
	reg check = 0;
	
	always @ (posedge ssend, negedge rst) begin
		if(~rst) begin
			sdata <= 0;
			check <= 0;
		end
		else begin
			sdata <= data;
			check <= ~check;
		end
	end
	
	always @ (posedge clk, negedge rst) begin
		if(~rst) begin
			i <= 0;
			dout = 1;
			state <= 2'h2;
			last_check = 0;
		end
		else begin
			if(check != last_check) begin
					state <= 2'h0;
					last_check = check;
			end
		
			case(state)
				2'h0: begin
					send <= 1;
					dout = 0;
					state <= 2'h1;
				end

				2'h1: begin
					dout = sdata[i];
					
					if(i != 7) begin
						i <= i + 1;
					end
					else begin
						i <= 0;
						state <= 2'h2;
					end
				end
					
				2'h2: begin
					dout = 1;
					send <= 0;
				end
					
			endcase
		end
	end
	
endmodule
Posted Updated Codes2 minutes read (About 309 words)

Verilog分频器与数码管

暑期课又在做FPGA了。

每次都写一遍分频器真的挺费事的，所以干脆现在给它记下来。

分频器

可以用于奇分频和偶分频。

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
module mydivider(
	input clk,
	input reset,
	output dclk);
	
	parameter width = 3;
	parameter n = 5;
	
	reg dclk_p, dclk_n;
	
	reg [width-1:0] count_p = 0;
	reg [width-1:0] count_n = 0;
	
	always @ (posedge clk, negedge reset) begin
		if(~reset) begin
			count_p <= 0;
			dclk_p <= 0;
		end
		else begin
			count_p <= count_p + 1;
			case(count_p)
				n/2-1: dclk_p <= 0;
				n-1: begin
					dclk_p <= 1;
					count_p <= 0;
				end
				default: dclk_p <= dclk_p;
			endcase
		end
	end
	
	always @ (negedge clk, negedge reset) begin
		if(~reset) begin
			count_n <= 0;
			dclk_n <= 0;
		end
		else begin
			count_n <= count_n + 1;
			case(count_n)
				n/2-1: dclk_n <= 0;
				n-1: begin
					dclk_n <= 1;
					count_n <= 0;
				end
				default: dclk_n <= dclk_n;
			endcase
		end
	end
	
	assign dclk = n%2 == 1 ? dclk_p | dclk_n : dclk_p;
	
endmodule

数码管 共阳

没有点的数码管代码（问就是板子没连，我也不知道为啥…

记录这个纯属是不想再写一遍这些管子的的亮灭了。

16个case分别对应0 - F的显示。

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
module nixietube(
	input [3:0] data,
	output reg [6:0] show);
	
	always begin
		case(data[3:0])
			4'd0: show <= 7'b1000000;
			4'd1: show <= 7'b1111001;
			4'd2: show <= 7'b0100100;
			4'd3: show <= 7'b0110000;
			4'd4: show <= 7'b0011001;
			4'd5: show <= 7'b0010010;
			4'd6: show <= 7'b0000010;
			4'd7: show <= 7'b1111000;
			4'd8: show <= 7'b0000000;
			4'd9: show <= 7'b0010000;
			4'd10: show <= 7'b0001000;
			4'd11: show <= 7'b0000011;
			4'd12: show <= 7'b1000110;
			4'd13: show <= 7'b0100001;
			4'd14: show <= 7'b0000110;
			4'd15: show <= 7'b0001110;
			default: show <= 7'b1000000;
		endcase
	end

endmodule

「向着星辰与深渊」

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 Codes11 minutes read (About 1633 words)

RTMP、HLS与HTTP-FLV

最近，有个项目需要使用OpenCV输出的图片推流到网页。这个东西相当的麻烦，主要是我以前从来没做过相关的东西。

这种需求相关的协议有三种：

协议 优点 缺点
RTSP 超低时延，可以达到毫秒级。 技术实现复杂。
RTMP 低时延，秒级；适应性好。 需要Flash支持，目前不流行。
HLS 动态切换码流；苹果安卓都可以用。 高时延，不适合做视频直播，一般用作点播或者音频广播。

首先排除RTSP，这东西几乎没啥可用文档。

  • RTMP

RTMP是应用层协议，全称Real Time Message Protocol。Adobe开发的协议，部署起来很简单，但是已经停止开发了。

协议细节请点击这里查看。

不过除非你要基于TCP/IP写一个RMTP应用，一般是不用太深入了解，只需要知道它能分组传输消息就行。

RTMP是基于FLV文件格式的协议，这种文件格式专门为流媒体设计，可以用一种持续不断“流”的方式进行传输。播放器拿到这种格式的分包后，会根据其中信息进行重排序，进而得到一个完整的媒体文件。

  • HLS

不同于RTMP，HLS同是应用层协议，但是其是基于HTTP协议传输的。

这个协议最大的特点就是把一个完整的媒体文件切分成许多短媒体，然后在一个m3u8文件里控制这些短媒体文件的读取。

本质上是在传输这些切分好的文件。

下面是CCTV1频道的语音主m3u8文件，它指出了各种带宽应使用的源。

依据奈氏准则，理想低通信道的最高码元传输速率 = 2W Baud

1
2
3
4
5
6
7
8
9
10
#EXTM3U
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=1800000,RESOLUTION=1280x720
/cctvwbcd/cdrmcctv1_1/index.m3u8?BR=td
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=1350000,RESOLUTION=1024x576
/cctvwbcd/cdrmcctv1_1/index.m3u8?BR=ud
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=900000,RESOLUTION=854x480
/cctvwbcd/cdrmcctv1_1/index.m3u8?BR=hd
#EXT-X-STREAM-INF:PROGRAM-ID=1,BANDWIDTH=600000,RESOLUTION=640x360
/cctvwbcd/cdrmcctv1_1/index.m3u8?BR=md

下面则是从主m3u8找到的对应的子m3u8，它控制着每个短流媒体的读取。

在直播场景下，子m3u8不能有#EXT-X-ENDLIST标志，否则浏览器会停止请求子m3u8，进而导致播放列表无法更新，直播断流。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#EXTM3U
#EXT-X-VERSION:3
#EXT-X-ALLOW-CACHE:NO
#EXT-X-TARGETDURATION:12
#EXT-X-MEDIA-SEQUENCE:1684822835
#EXTINF:12.000000,
#EXT-X-PROGRAM-DATE-TIME:2023-05-29T16:33:34+0800
cdrmcctv1_1_1800-1684822835.ts
#EXTINF:12.000000,
#EXT-X-PROGRAM-DATE-TIME:2023-05-29T16:33:46+0800
cdrmcctv1_1_1800-1684822836.ts
#EXTINF:12.000000,
#EXT-X-PROGRAM-DATE-TIME:2023-05-29T16:33:58+0800
cdrmcctv1_1_1800-1684822837.ts
#EXTINF:12.000000,
#EXT-X-PROGRAM-DATE-TIME:2023-05-29T16:34:10+0800
cdrmcctv1_1_1800-1684822838.ts

实践证明，HLS大多数时间都是用作点播，直播支持较差。

上面两个协议使用起来体验都不太好。

前者要求浏览器拥有Flash支持，几乎只能使用本地的VLC媒体播放器；后者则是要求对媒体进行切分，时间粒度较大，延时较高，而且直播支持不太完善。

所以，对于当下的直播，大多数使用的是HTTP-FLV。

  • HTTP-FLV

这种方式结合了RTMP和HLS的优点，使用HTTP来分发FLV流。

在浏览器的网络监视器中，可以看到一个持续不断的文件请求，类型是x-flv。

实际上，我们也可以用HTTPS来传输。

  • 服务器配置

Nginx有专门的包对RTMP和HTTP-FLV做了支持（HLS一般访问文件即可），所以使用起来很简单。

首先需要从GitHub上拉取nginx-http-flv-module。这个模块已经支持了RTMP，所以不用另外加上别的模块。

解压后执行下面的命令进行编译安装：

1
2
3
./configure --add-module=/tmp/nginx-http-flv-module --with-http_ssl_module
make
make install

安完后默认在/usr/local/nginx下，可执行文件在上述路径的sbin下。

我的网站有SSL证书，所以配置了HTTPS。参考配置如下：

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
worker_processes  1;

events {
    worker_connections  1024;
}
 
# 添加RTMP服务
rtmp {
    server {
        listen 1935;
        application live {  
            live on;
        }
    }
}

# 添加http-flv服务
http {
    include             mime.types;
    default_type        application/octet-stream;
    sendfile            on;
    keepalive_timeout   65;

    server {
        listen 443;
        server_name blueberrycat.site;

		ssl on;
		#ssl证书的pem文件路径
		ssl_certificate  /root/blueberrycat.site_nginx/blueberrycat.site_bundle.pem;
		#ssl证书的key文件路径
		ssl_certificate_key /root/blueberrycat.site_nginx/blueberrycat.site.key;

		location / {
            root html;
            index index.html index.htm;
        }

		error_page 500 502 503 504  /50x.html;
        location = /50x.html {
            root html;
        }
        
        location /live {
            flv_live on;
            chunked_transfer_encoding  on;
            # 添加一些控制访问的头
            add_header 'Access-Control-Allow-Credentials' 'true';
            add_header 'Access-Control-Allow-Origin' '*';
            add_header 'Access-Control-Allow-Headers' 'X-Requested-With';
            add_header 'Access-Control-Allow-Methods' 'GET,POST,OPTIONS';
            add_header 'Cache-Control' 'no-cache';
        }
    }

	server {
		listen 80;
		server_name blueberrycat.site;
		# 将请求转成https
		rewrite ^(.*)$ https://$host$1 permanent;
	}
}

  • 推流

推流地址可以是

1
https://blueberrycat.site/live?port=1935&app=live&stream={ 自定义的推流名 }

也可以是

1
rtmp://blueberrycat.site/live/{ 自定义的推流名 }

取决于你想用哪种协议。

回到项目，Python与OpenCV推流代码如下：

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
import cv2
import subprocess as sp

cap = cv2.VideoCapture(0)

fps = cap.get(cv2.CAP_PROP_FPS)

size = (int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)),
        int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT)))
sizeStr = str(size[0]) + 'x' + str(size[1])

command = ['ffmpeg',
           '-y',                                    # 覆盖确认
           '-f', 'rawvideo',                        # 指定输入格式
           '-pix_fmt', 'bgr24',                     # 格式
           '-s', sizeStr,                           # 大小
           '-r', str(fps),                          # FPS
           '-i', '-',                               # 输入-管道
           '-c:v', 'libx264',                       # 264编码
           '-pix_fmt', 'yuv420p',                   # 420p
           '-preset', 'ultrafast',                  # 高速处理
           '-f', 'flv',                             # flv输出
           '-b', '1000000',                         # 码率
           'rtmp://blueberrycat.site/live/123']     # 推流地址

pipe = sp.Popen(command, stdin=sp.PIPE, shell=False, bufsize=10**8)     # 建立管道，这步很重要

while cap.isOpened():
    ret, frame = cap.read()
    if ret == True:
        cv2.imshow('frame', frame)
        if cv2.waitKey(1) & 0xFF == ord('q'):
            break
    else:
        break

    pipe.stdin.write(frame.tobytes())       # 推流

为了测试这些协议，我在旧网站上做了大量的实验，结果直接把主机弄炸了，旧网站直接下线。

默哀。

别在 Ubuntu 上用 yum 装 g++

  • 网页端

建议使用hlv.js，Bilibili开发（难以置信）。

官方说明

An HTML5 Flash Video (FLV) Player written in pure JavaScript without Flash. LONG LIVE FLV!

😂

使用示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
<script src="flv.min.js"></script>
<video id="videoElement"></video>
<script>
    if (flvjs.isSupported()) {
        var videoElement = document.getElementById('videoElement');
        var flvPlayer = flvjs.createPlayer({
            type: 'flv',
            url: '{ 这里换成自己的推流 }'
        });
        flvPlayer.attachMediaElement(videoElement);
        flvPlayer.load();
        flvPlayer.play();
    }
</script>

最终效果是这样的

直播效果图

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
×