写出matlab的程序，对 BPSK、QPSK、8-PSK、16-QAM、16-PSK，计算当 Eb/N0（dB）从-1.59dB 以 0.2dB 为 间隔变化到 16dB 时，相应的可达信息速率（即所有输入符号等概使用），并作图。

由于题目中要求对不同的调制方式进行计算和作图，因此需要分别编写针对不同调制方式的程序。以下是分别针对BPSK、QPSK、8-PSK、16-QAM、16-PSK进行计算和作图的程序：

BPSK：

EbN0dB = -1.59:0.2:16; % Eb/N0的范围
EbN0 = 10.^(EbN0dB/10); % dB转换为倍数

% 计算每个Eb/N0下的误码率
BER = qfunc(sqrt(2*EbN0));

% 计算相应的可达信息速率
R = 1 - BER;

% 作图
semilogy(EbN0dB, R);
xlabel('Eb/N0 (dB)');
ylabel('可达信息速率');
title('BPSK');

QPSK：

EbN0dB = -1.59:0.2:16; % Eb/N0的范围
EbN0 = 10.^(EbN0dB/10); % dB转换为倍数

% 计算每个Eb/N0下的误码率
BER = qfunc(sqrt(EbN0));

% 计算相应的可达信息速率
R = 2*(1 - BER);

% 作图
semilogy(EbN0dB, R);
xlabel('Eb/N0 (dB)');
ylabel('可达信息速率');
title('QPSK');

8-PSK：

EbN0dB = -1.59:0.2:16; % Eb/N0的范围
EbN0 = 10.^(EbN0dB/10); % dB转换为倍数

% 计算每个Eb/N0下的误码率
BER = qfunc(sqrt(2*sin(pi/8)^2*EbN0));

% 计算相应的可达信息速率
R = 3*(1 - BER)*log2(8)/2;

% 作图
semilogy(EbN0dB, R);
xlabel('Eb/N0 (dB)');
ylabel('可达信息速率');
title('8-PSK');

16-QAM：

EbN0dB = -1.59:0.2:16; % Eb/N0的范围
EbN0 = 10.^(EbN0dB/10); % dB转换为倍数

% 计算每个Eb/N0下的误码率
BER = 3/2*qfunc(sqrt(4/10*EbN0));

% 计算相应的可达信息速率
R = 4*(1 - BER)*log2(16)/2;

% 作图
semilogy(EbN0dB, R);
xlabel('Eb/N0 (dB)');
ylabel('可达信息速率');
title('16-QAM');

16-PSK：

EbN0dB = -1.59:0.2:16; % Eb/N0的范围
EbN0 = 10.^(EbN0dB/10); % dB转换为倍数

% 计算每个Eb/N0下的误码率
BER = qfunc(sqrt(sin(pi/16)^2*EbN0));

% 计算相应的可达信息速率
R = 4*(1 - BER)*log2(16)/2;

% 作图
semilogy(EbN0dB, R);
xlabel('Eb/N0 (dB)');
ylabel('可达信息速率');
title('16-PSK');

需要注意的是，以上程序中用到了qfunc函数，这是MATLAB中的一个内置函数，用于计算Q函数的值。在计算QPSK和16-PSK的误码率时，用到了sin(pi/8)和sin(pi/16)这两个值，需要注意在MATLAB中这些值要用弧度计算。最后，程序通过调用semilogy函数绘制了以dB为横坐标、可达信息速率为纵坐标的图像。