matlab基本通信模块的建模(笔记)

求滤波器的最小阶和3dB截止频率

Matlab 中提供了 butter，cheb1ord，cheb2ord，ellipord 四个函数来分别设计巴特沃斯型、切比雪夫 1、2 型滤波器以及椭圆型模拟滤波器或数字滤波器。它们的调用格式相同，如下：

[n,Wn] = buttord(Wp,Ws,Rp,Rs) % 巴特沃斯型数字滤波器
[n,fn] = buttord(fp,fs,Rp,Rs,’s’) % 巴特沃斯型模拟滤波器
[n,Wn] = cheb1ord(Wp,Ws,Rp,Rs) % 切比雪夫1型数字滤波器
[n,fn] = cheb1ord(fp,fs,Rp,Rs,’s’) % 切比雪夫1型模拟滤波器
[n,Wn] = cheb2ord(Wp,Ws,Rp,Rs) % 切比雪夫2型数字滤波器
[n,fn] = cheb2ord(fp,fs,Rp,Rs,’s’) % 切比雪夫2型模拟滤波器
[n,Wn] = ellipord(Wp,Ws,Rp,Rs) % 椭圆型数字滤波器
[n,fn] = ellipord(fp,fs,Rp,Rs,’s’) % 椭圆型模拟滤波器

对于数字滤波器设计，输入参数 Wp，Ws 分别为归一化的频率 w p 和 w s 。对于模拟滤波器设计，输入参数 fp，fs 是不归一化的，即 f p 和 f s 。Rp，Rs 即以分贝为单位的通带内波动 R p 和阻带内最小衰减 R s 。返回值 n 为达到设计指标的最低系统阶数，Wn 为数字滤波器的 3dB 归一化截止频率，fn 为模拟滤波器的 3dB 截止频率。

• 低通数字滤波器情况：Wp<Ws，通带为 0 到 Wp，阻带为 Ws 到 1。
• 低通模拟滤波器情况：fp<fs，通带为 0 到 fp，阻带为 fs 到无穷。
• 高通数字滤波器情况：Wp>Ws，阻带为 0 到 Ws，通带为 Wp 到 1。
• 高通模拟滤波器情况：fp>fs，阻带为 0 到 fs，通带为 fp 到无穷。
• 带通数字滤波器情况：Ws(1)<Wp(1)<Wp(2)<Ws(2)，阻带为 0 到 Ws(1) 以及 Ws(2)
到 1，通带为 Wp(1) 到 Wp(2)。
• 带通模拟滤波器情况：fs(1)<fp(1)<fp(2)<fs(2)，阻带为 0 到 fs(1) 以及 fs(2) 到无穷，通带为 fp(1) 到 fp(2)。
• 带阻数字滤波器情况：Wp(1)<Ws(1)<Ws(2)<Wp(2)，通带为 0 到 Wp(1) 以及 Wp(2)到 1，阻带为 Ws(1) 到 Ws(2)。
• 带阻模拟滤波器情况：fp(1)<fs(1)<fs(2)<fp(2)，通带为 0 到 fp(1) 以及 fp(2) 到无穷，阻带为 fs(1) 到 fs(2)。

计算滤波器的传递函数

求出滤波器的阶数以及 3dB 截止频率後，可用相应的 Matlab 函数计算出实现传递函数的分子分母系数来。巴特沃斯型滤波器是通带内最大平坦、带外单调下降型的，其计算命令是：

[b,a] = butter(n,Wn) % 计算数字低通或带通情况
[b,a] = butter(n,Wn,’ftype’) % 计算数字高通或带阻情况
[b,a] = butter(n,Wn,’s’) % 计算模拟低通或带通情况
[b,a] = butter(n,Wn,’ftype’,’s’)% 计算模拟高通或带阻情况

其中，对於数字滤波器，Wn 就是 3dB 归一化截止频率。对於模拟滤波器，Wn 则是未归一化的角频率（单位 rad/s），与 fn 的关系是 Wn=2pifn。当截止频率参数为 2 个元素的向量时，为计算带通或带阻滤波器，否则是计算高通或低通滤波器的。当ftype 为 high 时为计算高通，当 ftype 为 stop 时为计算带阻。对於数字滤波器而言，返回值 b，a 分别是传递函数 H(z) 的分子和分母多项式的系数矩阵。对於模拟滤波器则返回值 b，a 分别是传递函数 H(s) 的分子和分母多项式的系数矩阵。
切比雪夫 1 型滤波器是通带等波纹（Equiripple）、阻带单调下降型的，其计算命令是：

[b,a] = cheby1(n,Rp,Wn) % 计算数字低通或带通情况
[b,a] = cheby1(n,Rp,Wn,’ftype’) % 计算数字高通或带阻情况
[b,a] = cheby1(n,Rp,Wn,’s’) % 计算模拟低通或带通情况
[b,a] = cheby1(n,Rp,Wn,’ftype’,’s’) % 计算模拟高通或带阻情况

其中，Rp 为带内波动参数（dB）。其余参数含义与巴特沃斯型滤波器计算函数中的参数相同。
切比雪夫 2 型滤波器是通带内单调、阻带等波纹型的，其计算命令是：

[b,a] = cheby2(n,Rs,Wn) % 计算数字低通或带通情况
[b,a] = cheby2(n,Rs,Wn,’ftype’) % 计算数字高通或带阻情况
[b,a] = cheby2(n,Rs,Wn,’s’) % 计算模拟低通或带通情况
[b,a] = cheby2(n,Rs,Wn,’ftype’,’s’) % 计算模拟高通或带阻情况

其中，Rs 为阻带内最小衰减参数（dB）。
椭圆型滤波器是通带、阻带内均为等波纹型的。其计算命令是：

[b,a] = ellip(n,Rp,Rs,Wn) % 计算数字低通或带通情况
[b,a] = ellip(n,Rp,Rs,Wn,’ftype’) % 计算数字高通或带阻情况
[b,a] = ellip(n,Rp,Rs,Wn,’s’) % 计算模拟低通或带通情况
[b,a] = ellip(n,Rp,Rs,Wn,’ftype’,’s’) % 计算模拟高通或带阻情况

其中，Rp 为通带内波动参数（dB），Rs 为阻带内最小衰减参数（dB）。

实例

例题1 设计一个模拟低通滤波器

试设计一个模拟低通滤波器，f p = 2400Hz，f s = 5000Hz，R p = 3 dB，R s =25dB。分别用巴特沃斯和椭圆滤波器原型，求出其 3dB 截止频率和滤波器阶数，传递函数，并作出幅频、相频特性曲线。

巴特沃斯滤波器设计的代码如下：

clear;
f_p=2400; f_s=5000; R_p=3; R_s=25;     % 设计要求指标
[n, fn]=buttord(f_p,f_s,R_p,R_s, 's'); % 计算阶数和截止频率
Wn=2*pi*fn;                            % 转换为角频率
[b,a]=butter(n, Wn, 's');              % 计算H(s)
f=0:100:10000;                         % 计算频率点和频率范围
s=j*2*pi*f;                            % s=jw=j*2*pi*f
H_s=polyval(b,s)./polyval(a,s);        % 计算相应频率点处H(s)的值
figure(1);
subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性
axis([0 10000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_s));         % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');
figure(2); freqs(b,a);  % 也可用指令freqs直接画出H(s)的频率响应曲线。

椭圆滤波器设计的程序代码如下：

% ch3example1B.m
clear;
f_p=2400; f_s=5000; R_p=3; R_s=25;     % 设计要求指标
[n, fn]=ellipord(f_p,f_s,R_p,R_s,'s'); % 计算阶数和截止频率
Wn=2*pi*fn;                            % 转换为角频率
[b,a]=ellip(n,R_p,R_s,Wn,'s');         % 计算H(s)
f=0:100:10000;                         % 计算频率点和频率范围
s=j*2*pi*f;                            % s=jw=j*2*pi*f
H_s=polyval(b,s)./polyval(a,s);        % 计算相应频率点处H(s)的值
figure(1);
subplot(2,1,1); plot(f, 20*log10(abs(H_s))); % 幅频特性
axis([0 10000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_s));         % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');
figure(2); freqs(b,a);   % 也可用指令freqs直接画出H(s)的频率响应曲线。

例题2 设计巴特沃斯数字低通滤波器

试设计一个巴特沃斯型数字低通滤波器，设采样率为 8000Hz，f p =2100Hz，f s = 2500Hz，R p = 3dB，R s = 25dB。

% ch3example2A.m
f_N=8000;                          % 采样率
f_p=2100; f_s=2500; R_p=3; R_s=25; % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2);    % 计算归一化频率
[n, Wn]=buttord(Wp,Ws,R_p,R_s);    % 计算阶数和截止频率
[b,a]=butter(n, Wn);               % 计算H(z)
figure(1);
freqz(b,a, 1000, 8000)             % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 4000 -30 3])
figure(2);                              % 第二种作图方法
f=0:40:4000;                            % 计算频率点和频率范围
z=exp(j*2*pi*f./(f_N));                 % 
H_z=polyval(b,z)./polyval(a,z);         % 计算相应频率点处H(s)的值
subplot(2,1,1); plot(f, 20*log10(abs(H_z))); % 幅频特性
axis([0 4000 -40 1]);
xlabel('频率 Hz');ylabel('幅度 dB');
subplot(2,1,2); plot(f, angle(H_z));         % 相频特性
xlabel('频率 Hz');ylabel('相角 rad');

例题3 切比雪夫数字高通滤波器

试设计一个切比雪夫 1 型高通数字滤波器，采样率为 8000Hz，f p =1000Hz，f s = 700Hz，R p = 3dB，R s = 20dB。

% ch3example3A.m
f_N=8000;                            % 采样率
f_p=1000; f_s=700; R_p=3; R_s=20;    % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2);      % 计算归一化频率
[n, Wn]=cheb1ord(Wp,Ws,R_p,R_s);     % 计算阶数和截止频率
[b,a]=cheby1(n, R_p, Wn, 'high');    % 计算H(z)
freqz(b,a, 1000, 8000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 4000 -30 3])

例题4 椭圆带通数字滤波器

试设计一椭圆型带通数字滤波器。设采样率为 10000Hz，f p = [1000,1500]Hz，f s = [600,1900]Hz，R p = 3dB，R s = 20dB。

% ch3example4A.m
f_N=10000;                                        % 采样率
f_p=[1000, 1500]; f_s=[600, 1900]; R_p=3; R_s=20; % 设计要求指标
Ws=f_s/(f_N/2); Wp=f_p/(f_N/2);                   % 计算归一化频率
[n, Wn]=ellipord(Wp,Ws,R_p,R_s);                  % 计算阶数和截止频率
[b,a]=ellip(n, R_p, R_s, Wn);                     % 计算H(z)
freqz(b,a, 1000, 10000) % 作出H(z)的幅频相频图, freqz(b,a, 计算点数, 采样率)
subplot(2,1,1); axis([0 5000 -30 3])

信源模型

例题1 求信号的近似频谱

求信号 f(t) = 2sin(2π100t) + cos(2π180t) 的近似频谱，要求计算频谱范围为 0 到 500Hz，频率分辨率（即频率点间隔）∆f 6 5Hz。用程序实现。

% ch3example12prg1.m
Df=5;      % 频率间隔
f_s=2*500; % 采样率
N=f_s/Df;  % 序列点数
t=0:1/f_s:(N-1)/f_s;    % 计算时间段
freq=0:Df:(N-1)*Df;       % 计算频率段
f_t=2*sin(2*pi*100*t)+cos(2*pi*180*t); % 信号
F_f=1/f_s*fft(f_t,N);     % 用FFT计算频谱
plot(freq-f_s/2,abs(fftshift(F_f)));   % 将零频率移动到FFT中心
xlabel('频率 Hz'); ylabel('幅度谱');   % 并作出幅度频谱

因为频谱的范围为0~500hz，所以连续信号的采样率为2×500hz。
使用fftshift()将计算频谱的中心位移到零频率，即输出频率的范围为[-500,500]Hz。

例题2 周期图法功率谱估计和计算信号功率

已知随机信号为
x(t) = sin2π50t + 2sin2π130t + n(t) (3.47)其中，n(t) 为零均值方差为 1 的高斯噪声。试用周期图法估计其功率谱密度。要求估计的频率范围为 [0,250]Hz，估计的频率分辨率为 1Hz。

% ch3example16prg1.m
fs=500;            % 采样率
Df=1;              % 频率分辨率
N=floor(fs/Df)+1;  % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs;         % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
                   % 截取时间段上的离散信号样点序列
Pxx=abs(fft(xk)).^2/(N^2);   % 功率谱估计
Pav_timedomaim=sum(xk.^2)/N  % 在时域计算信号功率
Pav_freqdomain=sum(Pxx)      % 通过功率谱计算信号功率
plot(F,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB') % 作出功率谱密度图

Pav_timedomaim =

    3.5432

Pav_freqdomain =

    3.5432

随着采样点数增加，该估计是渐近无偏的。从图中可以看出，采用周期图法估计得出的功率谱很不平滑，相应的估计协方差比较大。而且采用增加采样点的办法也不能使得周期图变得更加平滑，这是周期图法的缺点。
对周期图法进行改进的思想是将信号分段进行估计，然後再将这些估计结果进行平均，从而减小估计的协方差，使得估计功率谱图变得平滑。例如，将以上 501 点的信号分为 3 段，分别作周期图法估计，然後加以平均。程序代码如下。

% ch3example16prg2.m
fs=500;            % 采样率
Df=1;              % 频率分辨率
N=floor(fs/Df)+1;  % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs;         % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
                   % 截取时间段上的离散信号样点序列
Pxx=(abs(fft(xk(1:167))).^2+...
     abs(fft(xk(168:334))).^2+...
     abs(fft(xk(335:501))).^2)/3/((N/3)^2);
Pav_timedomaim=sum(xk.^2)/N  % 在时域计算信号功率
Pav_freqdomain=sum(Pxx)      % 通过功率谱计算信号功率
plot(0:3:fs,10*log10(Pxx));xlabel('freq Hz');ylabel('PSD dB');
%作出功率谱密度图

Matlab 中提供了多种估计功率谱的算法。最常用的是用指令「psd」、「pwelch」求功率谱。它们都使用 Welch 平均修正周期图法来进行谱估计。但这两个指令的使用方法稍有不同。指令「 psd 」的调用格式是：

Pxx= psd(X,NFFT,Fs,WINDOW,NOVERLAP);
[Pxx,F] = psd(X,NFFT,Fs,WINDOW,NOVERLAP);
其中:
pxx是功率谱密度(纵轴);
F为频率(横轴),F的点数是NFFT的一半;
X是样本信号序列;
NFFT是进行FFT的点数(默认为256);
Fs为指定的采样率(默认为2);
WINDOW是窗函数(默认为海明窗256点);
NOVERLAP是分段混叠的点数(默认为0)。

需要注意的是，「psd」函数的计算结果没有以 FFT 的点数进行归一化。所以进行归一化处理时需要再除以 FFT 点数的一半，也即，由功率谱密度结果序列 Pxx 来计算信号功率的方法是：
信号功率=sum(Pxx)./(FFT点数的一半)

例如，对信号序列 xk 采用 512 点 FFT，采样率为 500Hz，并使用海明窗 256 点，分段混叠点数为 128 点的 Welch 平均修正周期图法估计频谱的计算程序为：

% ch3example16prg5.m
fs=500;            % 采样率
Df=1;              % 频率分辨率
N=floor(fs/Df)+1;  % 计算的序列点数
t=0:1/fs:(N-1)/fs; % 截取信号的时间段
F=0:Df:fs;         % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
                   % 截取时间段上的离散信号样点序列
[Pxx,F] = psd(xk,512,500,hamming(256),128); 
plot(F,10*log10(Pxx/(512/2))); xlabel('freq Hz');ylabel('PSD dB')
Pav=sum(Pxx/(512/2)) % 通过功率谱计算信号功率
f:fs;         % 功率谱估计的频率分辨率和范围
xk=sin(2*pi*50*t)+2*sin(2*pi*130*t)+randn(1,length(t));
                   % 截取时间段上的离散信号样点序列
[Pxx,F] = psd(xk,512,500,hamming(256),128); 
plot(F,10*log10(Pxx/(512/2))); xlabel('freq Hz');ylabel('PSD dB')
Pav=sum(Pxx/(512/2)) % 通过功率谱计算信号功率

代码摘自书本《MATLAB/Simulink通信系统建模与仿真实例分析》