数字滤波器设计与信号处理技术

table {
border-collapse: collapse;
width: 100%;
margin-bottom: 1rem;
}
th, td {
border: 1px solid #ddd;
padding: 8px;
text-align: left;
}
th {
background-color: #f2f2f2;
}
tr:nth-child(even) {
background-color: #f9f9f9;
}
pre {
background-color: #f8f8f8;
padding: 15px;
border-radius: 4px;
overflow-x: auto;
}

65、设计一个低通滤波器，其规格如下：采样频率为2 Hz，通带上边缘频率为0.28 Hz，阻带低边缘频率为0.32 Hz，通带最大波纹为0.1 dB，阻带最小衰减为30 dB。

首先将相关参数转换为所需形式，如：

atmax = 0.1

（通带最大波纹）

atmin = 30

（阻带最小衰减）

Wp = 0.28 * 2 * pi

（通带上边缘频率，单位rad）

Ws = 0.32 * 2 * pi

（阻带低边缘频率，单位rad）

然后进行Butterworth模拟滤波器设计，计算阶数N：

N = ceil(log10((10^(0.1*atmin)-1)/(10^(0.1*atmax)-1))/Wnorm)

其中：

Wnorm = 2*log10(Ws/Wp)
Whp = Wp/((10^(.1*atmax)-1)^(1/(2*N)))
Omega = Whp

接着使用

buttap(N)

得到零点、极点和增益

[zc, pc, kc]

。最后进行频率变换：

pc = pc * Omega

不过，由于采样频率为2 Hz，在实际设计中可能还需要考虑数字滤波器设计的相关转换，可使用双线性变换等方法将模拟滤波器转换为数字滤波器。

66、设计一个低通滤波器，其规格如下：采样频率2Hz，通带上边缘0.28Hz，阻带 下边缘0.32Hz，通带最大波纹0.1dB，阻带最小衰减30dB。a. 找出满足这些规格的巴特沃斯、切比雪夫I型和II型以及椭圆滤波器的阶数。b. 为什么在相同的要求规格下，椭圆滤波器的阶数总是最低的？

a. 一般可使用以下MATLAB公式计算：

巴特沃斯

：

matlab
N = ceil(log10((10^(0.1*atmin)-1)/(10^(0.1*atmax)-1))/(2*log10(Ws/Wp)));


切比雪夫I型和II型、椭圆滤波器

：

可使用MATLAB函数

cheb1ord

、

cheb2ord

、

ellipord

计算。

b. 椭圆滤波器在通带和阻带都允许有波纹，它能以最少的阶数在给定频率范围内实现规定的衰减和波纹要求，因此在相同规格下，其阶数通常是最低的。

67、使用MATLAB函数boxcar、bartlett、hanning、hamming、blackman和kaiser绘制矩形、三角、汉宁、汉明、布莱克曼和凯泽窗。计算并绘制它们的频谱幅度。强调这些频谱的旁瓣电平与带宽。总结这些窗函数对有限长单位冲激响应（FIR）滤波器设计的影响。

不同窗函数在FIR滤波器设计中的影响

在FIR滤波器设计中，不同窗函数会对滤波器的性能产生不同影响：

矩形窗

：

主瓣最窄

频谱分辨率最佳

但旁瓣最高

动态分辨率最差

其他窗函数

（如汉宁窗、汉明窗、布莱克曼窗和凯泽窗）：

可改善动态分辨率

但会降低频谱分辨率

在设计FIR滤波器时，需要在以下因素之间进行权衡：

旁瓣电平

主瓣宽度

旁瓣衰减率

68、使用最小二乘法（LS 方法）计算一个 25 阶导数滤波器的系数，采样频率 fs = 20 kHz。绘制其传递函数、脉冲响应和零点。设想该滤波器的一个可能应用。

需按以下步骤求解：

首先根据最小二乘法原理和 25 阶导数滤波器、采样频率 $ f_s = 20 , ext{kHz} $ 的条件计算滤波器系数；

然后使用相关工具（如 MATLAB）绘制传递函数、脉冲响应和零点；

可能的应用例如在信号处理中检测信号的变化率，如在音频处理中检测声音信号的突变等。

69、使用频率采样法设计一个120阶的FIR高通滤波器。考虑截止频率fc = 6 kHz，采样频率fs = 20 kHz。绘制滤波器的传递函数。将其阶数与具有类似传递函数的IIR滤波器的阶数进行比较。

设计该滤波器需按频率采样法步骤确定离散频率点的期望频率响应，再用 IDFT 得到滤波器系数。绘制传递函数可借助 MATLAB 的

freqz

函数。

通常，达到相同滤波效果时，IIR 滤波器阶数低于 FIR 滤波器，因为 IIR 滤波器存在反馈环节，可利用极点实现更陡峭的过渡带。

70、a. 生成以下三种信号，每个信号包含1024个样本： – 两个归一化频率分别为0.075和0.175的正弦波之和； – 长度为50点的脉冲信号； – 由方差为1的零均值白高斯噪声驱动的带通椭圆滤波器的输出信号，滤波器规格如下：通带频率边缘为0.205和0.245，阻带频率边缘为0.195和0.255，通带波纹为0.5 dB，阻带最小衰减为50 dB。 向每个信号添加方差为3的零均值白高斯噪声。 b. 使用相应的维纳滤波器对三个含噪信号进行滤波。解释每种情况下的滤波效果。

a. 信号生成代码如下：

nivbr=3;
s1=sin(2*pi*.075*[1:1024])+sin(2*pi*.175*[1:1024]);
sf1=abs(fft(s1));
s2=[ones(1,50) zeros(1,1024-50)];
sf2=abs(fft(s2));
[n,wn] = ellipord([0.205,0.245],[0.195,0.255],0.5,50);
[b0,a0] = ellip(n,0.5,50,wn);
[h0,w0]=freqz(b0,a0,512);
s3=filter(b0,a0,randn(1,1024));
sf3=abs(fft(s3));
% 噪声添加
x1=s1+nivbr*randn(1,1024);
x1f=abs(fft(x1));
x2=s2+nivbr*randn(1,1024);
x2f=abs(fft(x2));
x3=s3+nivbr*randn(1,1024);
x3f=abs(fft(x3));

b. 维纳滤波代码如下：

nb=9;
na=10;
% 滤波器分子和分母阶数
hz=fft(xcorr(x1,s1))./fft(xcorr(x1));
[b,a]=invfreqz(hz,[1:2047]*2*pi/2047,nb,na);
[h1,w1]=freqz(b,a,512);
sr1=filter(b,a,x1);
srf1=abs(fft(sr1));
hz=fft(xcorr(x2,s2))./fft(xcorr(x2));
[b,a]=invfreqz(hz,[1:2047]*2*pi/2047,nb,na);
[h2,w2]=freqz(b,a,512);
sr2=filter(b,a,x2);
srf2=abs(fft(sr2));
hz=fft(xcorr(x3,s3))./fft(xcorr(x3));
[b,a]=invfreqz(hz,[1:2047]*2*pi/2047,nb,na);
[h3,w3]=freqz(b,a,512);
sr3=filter(b,a,x3);
srf3=abs(fft(sr3));

71、使用不同的基本脉冲（非归零双极性码、三角脉冲、曼彻斯特码）对发射波形进行信号生成，并重复类似的操作。比较大量仿真得到的理论和实验错误概率。

可按照以下步骤解决：

分别使用非归零双极性码、三角脉冲和曼彻斯特码作为基本脉冲生成信号；

对每种脉冲类型进行大量仿真，计算理论和实验错误概率；

比较不同脉冲类型下的理论和实验错误概率。

72、一个有源早期预警雷达以振幅A发射脉冲串（每秒1200个脉冲），使用均匀旋转天线（每分钟15转），主瓣宽度为1.5°。1. 被扫描目标会反向散射多少个脉冲？2. 通过适当的处理，被扫描目标反向散射的N个脉冲形成一个具有N个分量的观测向量，该向量受到标准差为σ的零均值加性高斯白噪声的干扰。使用奈曼 – 皮尔逊准则进行检测，该准则固定虚警概率水平并最大化检测概率。证明该准则等效于将N个分量的和与一个应定义的阈值进行比较。计算当A = 1，σ = 0.6，N = 20，虚警概率PFA = 10^(-10)时的结果。编写一个MATLAB代码来预测不同参数值下的检测性能，并绘制ROC（接收机工作特性）曲线。这些曲线展示了不同信噪比下检测概率随虚警概率的变化情况。

要解决此问题，可按以下步骤进行：

计算被扫描目标反向散射的脉冲数

：

先算出天线旋转一圈所需时间，再根据主瓣宽度算出主瓣扫描目标的时间，最后结合脉冲发射速率得出反向散射的脉冲数。

证明奈曼 – 皮尔逊准则等效于将N个分量的和与阈值比较

：

根据准则的原理和观测向量的特性进行推导。

计算给定参数下的结果

：

代入参数进行计算。

编写MATLAB代码

：

根据上述原理编写代码来预测不同参数下的检测性能并绘制ROC曲线。

73、给定一个高斯白过程，计算并绘制其平均周期图和修正周期图。找出它们的偏差和方差。改变平均序列的数量，并使用 MATLAB 函数 boxcar、triang、hamming、hanning 和 blackman 研究不同窗函数对修正周期图的影响。

可按以下步骤操作：

计算平均周期图

：使用无序列恢复的 MATLAB 代码进行计算；

计算修正周期图

：在进行 FFT 计算前，将每个序列乘以加权窗函数，可使用如

blackman

等窗函数；

找出偏差和方差

：对于平均周期图，该估计仍有偏差，但方差会随平均序列数量增加而减小；

改变平均序列数量并研究不同窗函数影响

：使用 MATLAB 函数

boxcar

、

triang

、

hamming

、

hanning

和

blackman

改变窗函数，观察对修正周期图的影响。

74、已知一个信号，使用修正周期图法对该信号进行频谱分析。先使用汉明窗，再使用布莱克曼窗，其中K = 512，L = 4。将此结果与之前获得的结果进行比较。使用修正周期图的算法计算信号的修正周期图。

以下是使用汉明窗和布莱克曼窗计算信号修正周期图的代码示例：

汉明窗

：

K = 128;
weighting_window = [hamming(K)];
weighting_window = weighting_window * sqrt(K / sum(weighting_window .* weighting_window));
weighting_window = weighting_window.';
[PSD_modified_periodogram, frequency] = modified_periodogram(signal, K, weighting_window);

布莱克曼窗

：

figure;
semilogy(frequency, PSD_modified_periodogram, 'b');
hold on;
weighting_window = [blackman(K)];
weighting_window = weighting_window * sqrt(K / sum(weighting_window .* weighting_window));
weighting_window = weighting_window.';
[PSD_modified_periodogram, frequency] = modified_periodogram(signal, K, weighting_window);
semilogy(frequency, PSD_modified_periodogram, 'r:');
xlabel('normalized frequency');
ylabel('Amplitude');
title('Signal modified periodogram calculated for N=512 and L=4');
legend('Hamming window', 'Blackman window');

与平均周期图相比，修正周期图中对应两个频谱分量的峰值更大，这是由于加权窗的主瓣宽度。这对于检测峰值的频率估计可能是一个缺点，但当这些频率不是谐波时可能有用。

75、生成以下两个频率非谐波关系的正弦波的混合信号，并添加方差为 0.031 的零均值白高斯噪声。第一个正弦波：频率 = 25.4 Hz，振幅 = 1；第二个正弦波：频率 = 51.3 Hz，振幅 = 0.01。计算该信号在信号长度 K 分别为 256、512、1024 时的周期图，并确定作为 K 的函数的谐波频率。然后对 {K = 2048, L = 8}、{K = 4096, L = 16} 和 {K = 4096, L = 8} 这几种情况测试平均周期图。在每种情况下，谐波频率分别是多少？它们是否依赖于 K？对这些结果进行评论并突出频谱泄漏现象。最后，对不同窗口应用修改后的周期图并得出结论。

需按要求编写代码进行信号生成、周期图计算、平均周期图测试、频谱泄漏分析以及修改后周期图应用等操作来得出结果。

76、给定一个信号，使用参数谱分析方法（假设信号为自回归（AR）模型）绘制其功率谱密度，并确定AR模型的阶数。将估计的功率谱密度与使用平均周期图法和修正周期图法得到的结果进行比较。

确定AR模型阶数，用参数谱分析绘制功率谱密度，获取平均和修正周期图法的结果进行比较。

77、生成以下两个频率非谐波关系的正弦波的混合信号，并添加方差为 0.031 的零均值白高斯噪声。第一个正弦波：频率 = 25.4 Hz，振幅 = 1；第二个正弦波：频率 = 51.3 Hz，振幅 = 0.01。计算该信号在信号长度 K 分别为 256、512、1024 时的周期图，并确定谐波频率与 K 的函数关系。然后对 {K = 2048, L = 8}、{K = 4096, L = 16} 和 {K = 4096, L = 8} 这几种情况测试平均周期图。每种情况下的谐波频率是多少？它们是否依赖于 K？对这些结果进行评论，并指出频谱泄漏现象。最后，对不同窗口应用修正周期图并得出结论。

需根据相关信号处理知识和算法进行计算和分析。一般步骤为：

先按要求生成混合信号并添加噪声；

再用相应算法计算不同 K 值下的周期图、特定情况的平均周期图和不同窗口的修正周期图；

分析谐波频率与 K 的关系、频谱泄漏现象等并得出结论。

78、使用 Haar 小波对含噪正弦信号进行多分辨率分析。然后验证合成互逆算法能够从最后一个近似级别和所有细节中重建信号。证明在重建之前通过消除部分分解系数可以对信号进行去噪。

要完成该任务，首先可使用

analyzehaar

函数对含噪正弦信号进行多分辨率分析得到近似和细节系数；接着实现合成互逆算法，利用最后近似级别和所有细节系数重建信号；最后在重建前消除部分分解系数，对比去噪前后信号，证明去噪可行性。

79、练习中提出计算wigner函数的第一步是计算与所考虑的实信号相关的解析信号。通过这种方式，频谱支撑被除以2，从而避免了混叠现象。通过将此过程与直接计算维格纳 – 威利分布（不使用解析信号）进行比较，来说明该过程的优势。考虑一个包含多个时频原子、以奈奎斯特频率采样的实信号。

可按以下思路解答：

先生成一个包含多个时频原子且以奈奎斯特频率采样的实信号；

然后分别用直接计算和先计算解析信号再计算的方法得到维格纳 – 威利分布；

最后对比两种方法的结果，分析使用解析信号避免混叠现象的优势。

80、编写一个 MATLAB 代码来计算与其他核（如 Choi – Williams、Rihaczek 等）相关的时频分布，并验证它们的理论特性。

编写此类代码可按以下步骤进行：

明确 Choi – Williams、Rihaczek 等核函数的数学定义；

生成测试信号；

根据核函数的定义编写计算时频分布的代码；

编写验证理论特性的代码，如验证时频聚集性、边缘特性等。

以下是一个简单的示例框架：

% 生成测试信号
N = 1024;
T = 1;
t = linspace(0, T, N);
x = sin(2*pi*10*t);

% Choi - Williams 时频分布计算示例（这里只是框架，需补充具体核函数计算）
function cw_distribution = choi_williams(x, t, N)
    cw_distribution = zeros(N, N);
    % 补充 Choi - Williams 核函数计算过程
    return;
end

cw = choi_williams(x, t, N);

% Rihaczek 时频分布计算示例（这里只是框架，需补充具体核函数计算）
function rihaczek_distribution = rihaczek(x, t, N)
    rihaczek_distribution = zeros(N, N);
    % 补充 Rihaczek 核函数计算过程
    return;
end

ri = rihaczek(x, t, N);

% 验证理论特性（这里只是框架，需补充具体验证方法）
function verify_properties(distribution)
    % 补充验证理论特性的代码
    return;
end

verify_properties(cw);
verify_properties(ri);

81、编写一个MATLAB函数，用于计算给定阶数的信号的分数阶傅里叶变换（FRFT）。

分数阶傅里叶变换（FRFT）MATLAB函数实现

以下是计算FRFT的MATLAB代码，已封装成函数：

function FRFT = FracFR(x0, a)
    deltax = sqrt(length(x0)/2);
    phi = a*pi/2;
    N = fix(length(x0)/2);
    deltax1 = 2*deltax;
    alpha = 1/tan(phi);
    beta = 1/sin(phi);
    % Kernel implementation
    T = [-N:N-1]/deltax1;
    T = T(:);
    x0 = x0(1:2*N);
    f1 = exp(-i*pi*tan(phi/2)*T.*T);
    f1 = f1(:);
    x = x0.*f1;
    clear T;
    t = [-2*N+1:2*N-1]/deltax1;
    hlptc = exp(i*pi*beta*t.*t);
    clear t;
    hlptc = hlptc(:);
    N2 = length(hlptc);
    N3 = 2^(ceil(log(N2+2*N-1)/log(2)));
    hlptcz = [hlptc;zeros(N3 - N2,1)];
    xz = [x;zeros(N3-2*N,1)];
    Hcfft = ifft(fft(xz).*fft(hlptcz));
    clear hlptcz;
    clear xz;
    Hc = Hcfft(2*N:4*N-1);
    clear Hcfft;
    clear hlptc;
    Aphi = exp(-i*(pi*sign(sin(phi))/4-phi/2))/sqrt(abs(sin(phi)));
    xx = [-N:N-1]/deltax1;
    f1 = f1(:);
    FRFT = (Aphi*f1.*Hc)/deltax1;
end

使用示例

% 示例信号
x0 = randn(100, 1);

% 阶数
a = 0.5;

% 计算FRFT
FRFT = FracFR(x0, a);

82、编写一个MATLAB函数，根据给定的输入信号计算该信号的M阶矩。

一般而言，在 MATLAB 中可以通过以下思路编写函数来计算 M 阶矩：

首先明确 M 阶矩的计算公式，对于离散信号，M 阶矩通常是信号值的 M 次方的平均值。

以下是一个简单示例代码：

function moment = mth_order_moment(signal, M)
    moment = mean(signal.^M);
end

此代码定义了一个名为

mth_order_moment

的函数，接收信号和阶数 M 作为输入，返回该信号的 M 阶矩。

83、考虑以下信号：(s(t) = exp[j2π(0.1t + 0.0012t^2)] + exp[j2π(0.12t + 0.0012t^2)], t = 0..255) a. 绘制信号每个分量的瞬时频率规律。b. 使用分数阶傅里叶变换（FRFT）估计每个分量的参数。c. 提出提高所得分辨率的解决方案。

a. 可先分别对信号的两个分量求瞬时频率公式，再用绘图工具绘制；

b. 可通过对信号应用不同阶数的FRFT，找到使频谱能量最集中的阶数来估计参数；

c. 提高分辨率可增加采样点数、减小阶数的步长、使用更合适的窗函数等。

84、考虑一个具有两个相似分量的信号：(s(t)=exp[j2π0.02t^{1.5}]+exp[j2π0.01t^{1.5}], t = 0..255)。使用扭曲原理估计每个分量的调制率。计算该信号对应的非线性匹配时频表示（TFR）。将所得结果与两个分量的瞬时频率定律进行比较。

可使用一组扭曲算子并选择使预扭曲信号对应频谱峰宽度最小的算子，其阶数 $ k $ 表征信号调制，以此估计调制率。

计算非线性匹配 TFR 可采用针对双曲调制的计算方法。

比较结果时先求出两个分量的瞬时频率定律，再与计算得到的 TFR 结果进行对比。

85、a. 生成两类各200个二维高斯向量，其均值向量和协方差矩阵如下：第一类均值向量为[5; 3]，协方差矩阵为[[4, 6]; [6, 12]]；第二类均值向量为[4; 15]，协方差矩阵为[[12, -4]; [-4, 8]]。然后将每一类随机划分为训练子集和测试子集，且两个子集包含相同数量的向量。b. 假设两类的先验概率相等。使用从训练子集得到的贝叶斯分类器的直接近似对测试子集进行分类。c. 使用模糊K近邻（KNN）方法重复分类过程。在同一图中绘制测试向量以及两个分类器生成的分离面。d. 通过对1000次结果的性能进行平均来评估两个分类器的分类率。

a. 可使用以下MATLAB代码生成向量并划分训练集和测试集：

class_nbr = 2;
vect_nbr = 200;
mean_vect1 = [5; 3];
cov_matr1 = [4 6; 6 12];
class1 = mvnrnd(mean_vect1', cov_matr1, vect_nbr)';
mean_vect2 = [4; 15];
cov_matr2 = [12 -4; -4 8];
class2 = mvnrnd(mean_vect2', cov_matr2, vect_nbr)';
idx = randperm(vect_nbr);
vect_nbr_test = fix(vect_nbr / 2);
idx_test = idx(1:vect_nbr_test);
idx_train = idx(vect_nbr_test + 1:end);
vect_label_train1 = ones(1, vect_nbr - vect_nbr_test);
vect_label_test1 = ones(1, vect_nbr_test);
vect_label_train2 = 2 * vect_label_train1;
vect_label_test2 = 2 * vect_label_test1;
vect_train1 = class1(:, idx_train);
vect_test1 = class1(:, idx_test);
vect_train2 = class2(:, idx_train);
vect_test2 = class2(:, idx_test);
matr_train = [vect_train1 vect_train2];
vect_label_train = [vect_label_train1 vect_label_train2];
matr_test = [vect_test1 vect_test2];

b. 未给出此步骤具体代码。

c. 未给出此步骤具体代码。

d. 贝叶斯分类器对第一类和第二类的平均分类率分别为 0.969 和 0.979；模糊 KNN 分类器对第一类和第二类的平均分类率分别为 0.967 和 0.973。

86、a. 生成两类向量，每类1500个。第一类向量均匀分布在半径为3的圆盘内，第二类向量均匀分布在半径从3到6的圆环区域内。这两类向量的中心都在原点。然后将每类向量随机分成三个子集（训练集、测试集和验证集），每个子集包含相同数量的向量。绘制这两类的训练向量。b. 使用第一个子集训练一个多层感知器，并使用第三个子集在其泛化能力达到最大时决定训练过程的结束。考虑带有动量项和可变学习率的反向传播学习算法。使用训练好的多层感知器对测试向量进行分类。c. 绘制多层感知器对应于两类向量所占区域的输出，并强调其与向量分布的联系。

a部分给出了生成向量并划分训练、测试和验证子集的MATLAB代码的提示；c部分给出了绘制多层感知器对应两类区域输出的MATLAB代码：

figure; 
imagesc(xv,yv,g1_pmc); 
hold on; 
contour(xv,yv,g1_pmc); 
axis xy; 
colormap(spring); 
xlabel('x_1'); 
ylabel('x_2'); 
colorbar; 
axis equal; 
title('MLP output corresponding to the class 1');

figure; 
imagesc(xv,yv,g2_pmc); 
hold on; 
contour(xv,yv,g2_pmc); 
axis xy; 
colormap(spring); 
xlabel('x_1'); 
ylabel('x_2'); 
colorbar; 
axis equal; 
title('MLP output corresponding to the class 2');