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;
}
19、如果 N 是质数且 ρ 是 N 的原根，函数 ρ^i mod N（其中 i 取值范围是 1 到 N – 1）会将序列 {1, 2, 3, … , N – 1} 排列成序列 {ρ^1 mod N, ρ^2 mod N, ρ^3 mod N, …, ρ^(N – 1) mod N}。编写一个 MATLAB 代码来执行这些排列。
以下是实现该排列的 MATLAB 代码：
N = 5; % 示例 N 值
rho = 2; % 示例 rho 值
% 生成原始序列
original_sequence = 1:N-1;
% 执行排列
permuted_sequence = mod(rho.^(original_sequence-1), N);
% 显示结果
fprintf('原始序列: ');
disp(original_sequence);
fprintf('排列后的序列: ');
disp(permuted_sequence);
在上述代码中，首先定义了
N
和
rho
的值。然后生成原始序列
original_sequence
，接着使用
mod
函数执行排列操作，最后显示原始序列和排列后的序列。你可以根据需要修改
N
和
rho
的值。
20、在拉普拉斯域中定义一个抗混叠滤波器的传递函数，该滤波器在频率3500 Hz处衰减0.5 dB，在频率4500 Hz处衰减30 dB。
Fp = 3500;
Fs = 4500;
Wp = 2*pi*Fp;
Ws = 2*pi*Fs;
[N, Wn] = buttord(Wp, Ws, 0.5, 30,'s');
[b,a] = butter(N, Wn, 's');
wa = 0:(3*Ws)/511:3*Ws;
h = freqs(b,a,wa);
plot(wa/(2*pi), 20*log10(abs(h)));
grid;
xlabel('Frequency [Hz]');
ylabel('Gain [dB]');
title('Frequency response');
axis([0 3*Fs -60 5]);
21、使用函数 filter.m 计算由以下方程定义的系统的脉冲响应 h(n)：y(n) – 1.7654y(n – 1)+ 0.81y(n – 2)= x(n)+ 0.5x(n – 1)。绘制 -10 到 60 之间的 h(n)。将得到的结果与理论结果进行比较。
以下是实现该功能的 MATLAB 代码：
a = [1 -1.7654 0.81];
b = [1 0.5];
delta = [1 zeros(1,99)];
h = filter(b, a, delta);
n = [-10:59];
h1 = [zeros(1,10) h];
stem(n, h1(1:70));
xlabel('n');
ylabel('h(n)');
grid
运行代码后可得到 -10 到 60 之间的 h(n) 图像，可将其与理论结果进行比较。
22、计算截止频率为 1 kHz 的 5 阶模拟巴特沃斯滤波器的极点。
% 原型滤波器设计
N = 5;
[zc, pc, kc] = buttap(N);
% 频率变换
Omega = 2 * pi * 10^3;
pc = pc * Omega;
% 绘制极点
zplane(zc, pc)
23、设计一个采样频率为 1Hz 的 IIR 低通滤波器，在频率范围 [0 0.1306]Hz 内衰减为 0.75dB，阻带衰减至少为 20dB，阻带下限截止频率为 0.201Hz。使用巴特沃斯模型计算相关模拟滤波器的极点。
% wp: 通带上限截止频率 [rad]
% ws: 阻带下限截止频率 [rad]
% atmax: 通带波纹或最大允许通带损耗 [dB]
% atmin: 最小阻带衰减 [dB]
atmax = 0.75;
atmin = 20;
Wp = 0.1306 * 2 * pi;
Ws = 0.201 * 2 * pi;
% 巴特沃斯模拟滤波器设计
Wnorm = 2 * log10(Ws / Wp);
N = ceil(log10((10^(0.1 * atmin) - 1) / (10^(0.1 * atmax) - 1)) / Wnorm);
Whp = Wp / ((10^(0.1 * atmax) - 1)^(1 / (2 * N)));
Omega = Whp;
[zc, pc, kc] = buttap(N);
% 频率变换
pc = pc * Omega;
% 绘制极点
zplane(zc, pc);
24、考虑一个6阶椭圆滤波器，其归一化截止频率为0.2，通带波纹为0.05 dB，阻带最小衰减为60 dB。编写一个MATLAB代码，以说明5位量化对滤波器特性（传递函数、极点和零点）的影响。
在MATLAB中，原答案里
a2dT
函数并不存在，我们可以使用
round
函数来实现5位量化。以下是修正后的代码：
% 设计椭圆滤波器
[b,a] = ellip(6,0.05,60,2*0.2);
% 定义频率范围
w = 0:pi/255:pi;
% 计算未量化滤波器的频率响应
h = freqz(b,a,w);
% 计算增益
g = 20*log10(abs(h));
% 对滤波器系数进行5位量化
bq = round(b * 2^5) / 2^5;
aq = round(a * 2^5) / 2^5;
% 计算量化后滤波器的频率响应
hq = freqz(bq,aq,w);
% 计算量化后的增益
gq = 20*log10(abs(hq));
% 绘制频率响应曲线
subplot(2,1,1);
plot(w/(2*pi),g,'b',w/(2*pi),gq,'r--');
axis([0 0.5 -80 10]);
grid on;
xlabel('Normalized frequency');
ylabel('Gain [dB]');
legend('Original', 'Quantized');
% 计算并显示极点和零点
z = roots(b);
p = roots(a);
[zq,pq] = zpkdata(tf(bq,aq),'v');
subplot(2,1,2);
zplane(z,p);
hold on;
zplane(zq,pq);
legend('Original Zeros and Poles', 'Quantized Zeros and Poles');
这段代码首先设计了一个6阶椭圆滤波器，然后计算了未量化和5位量化后的滤波器频率响应，并绘制了增益曲线。同时，还计算并绘制了未量化和量化后的极点和零点。
25、使用窗函数法设计一个25阶的FIR低通滤波器，使其过渡带最小。考虑采样频率fs = 20 kHz，通带频率上下限分别为f1 = 3 kHz和f2 = 7 kHz。绘制所设计滤波器的传递函数、脉冲响应和零点图。
N = 25;
f1 = 3e3;
f2 = 7e3;
fs = 2e4;
wn1 = 2*f1/fs;
wn2 = 2*f2/fs;
b = fir1(N, [wn1 wn2], boxcar(N+1));
[h, w] = freqz(b, 1);
fq = fs*w/(2*pi);
amp = 20*log10(abs(h));
phase = unwrap(angle(h))*180/pi;
figure;
subplot(221);
plot(fq, amp);
axis([0 max(fq) min(amp) 10]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
grid;
subplot(223);
plot(fq, phase);
axis([0 max(fq) 1.2*min(phase) 1.2*max(phase)]);
xlabel('Frequency [Hz]');
ylabel('Phase [°]');
grid;
subplot(2,2,2);
impz(b, 1);
xlabel('n');
subplot(2,2,4);
zplane(b, 1);
26、使用频率采样法设计一个24阶的FIR高通滤波器。考虑截止频率fc = 6 kHz，采样频率fs = 20 kHz。绘制所设计滤波器的传递函数、脉冲响应和零点。
以下是实现该设计的MATLAB代码：
N = 24;
fc = 6e3;
fs = 20e3;
inc = 0.01;
wcn = 2 * fc / fs;
ff = [0:inc:1];
hh = [zeros(1, wcn / inc - 1) ones(1, (1 - wcn) / inc + 2)];
b = fir2(N, ff, hh);
[h, w] = freqz(b, 1);
fq = fs * w / (2 * pi);
amp = 20 * log10(abs(h));
phase = unwrap(angle(h)) * 180 / pi;
subplot(221);
plot(fq, amp);
axis([0 max(fq) min(amp) 10]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
grid;
subplot(223);
plot(fq, phase);
axis([0 max(fq) 1.2 * min(phase) 1.2 * max(phase)]);
xlabel('Frequency [Hz]');
ylabel('Phase [°]');
grid;
subplot(2, 2, 2);
impz(b, 1);
xlabel('n');
subplot(2, 2, 4);
zplane(b, 1);
运行上述代码可得到滤波器的传递函数、脉冲响应和零点的绘图。
27、使用最小二乘法（LS方法）计算一个25阶有限长单位冲激响应（FIR）带阻滤波器的系数。考虑采样频率fs = 20 kHz，以及以下阻带频率的上下边界：f1 = 3 kHz和f2 = 7 kHz。绘制所设计滤波器的传递函数、冲激响应和零点。设想这种类型滤波器的一个可能应用。
代码实现如下：
N=25;
f1=3e3;
f2=7e3;
fs=2e4;
wn1=2*f1/fs;
wn2=2*f2/fs;
inc=0.01;
ff=[0:inc:1-inc];
hh=[ones(1,round((wn1/inc)-1)) zeros(1,round(((wn2-wn1)/inc)+1))];
hh=[hh ones(1,round((1-wn2)/inc))];
b=firls(N,ff,hh);
[h,w] = freqz(b,1);
fq = fs*w/(2*pi);
amp=20*log10(abs(h));
phase=unwrap(angle(h))*180/pi;
figure;
subplot(221);
plot(fq,amp);
axis([0 max(fq) min(amp) 10])
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]') ;
grid
subplot(223);
plot(fq,phase);
axis([0 max(fq) 1.2*min(phase) 1.2*max(phase)])
xlabel('Frequency [Hz]');
ylabel('Phase [°]');
grid
subplot(2,2,2);
impz(b,1);
xlabel('n')
subplot(2,2,4);
zplane(b,1);
此滤波器虽在阻带不是等波纹的，但可用于抑制阻带内干扰或扰动，避免其破坏有用信号分量。若要在阻带获得更好的抑制效果，可使用更高阶的FIR滤波器或无限长单位冲激响应（IIR）滤波器。
28、使用 Remez 方法设计一个带通 FIR 滤波器，其中滤波器的参数为：第一个阻带截止频率范围是 2.5kHz – 3.5kHz，第二个阻带截止频率范围是 6.5kHz – 7.5kHz，采样频率为 20kHz。绘制所得到滤波器的传递函数、脉冲响应和零点，并与其他设计方法得到的滤波器进行比较。
f11=2.5e3;
f12=3.5e3;
f21=6.5e3;
f22=7.5e3;
fs=2e4;
wvect=[f11 f12 f21 f22];
[n,fo,mo,w] = remezord(wvect,[0 1 0],[0.04 0.03 0.05],fs);
b = remez(n,fo,mo,w);
[h,w] = freqz(b,1);
fq = fs*w/(2*pi);
amp=20*log10(abs(h));
phase=unwrap(angle(h))*180/pi;
figure;
subplot(221);
plot(fq,amp);
axis([0 max(fq) min(amp) 10]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
grid;
subplot(223);
plot(fq,phase);
axis([0 max(fq) 1.2*min(phase) 1.2*max(phase)]);
xlabel('Frequency [Hz]');
ylabel('Phase [°]');
grid;
subplot(2,2,2);
impz(b,1);
xlabel('n');
subplot(2,2,4);
zplane(b,1);
与其他设计方法得到的滤波器相比，此滤波器阶数为 25，可从其脉冲响应的样本数或零点数获取。该滤波器具有更好的阻带抑制（衰减增益约为 10 dB）和更小的通带波纹，且用 Remez 方法设计的滤波器在阻带是等波纹的。
29、考虑两个频率分别为1kHz和1.56kHz的正弦波之和，采样频率为10kHz。使用Remez方法设计一个低通FIR滤波器，从该混合信号中提取第一个正弦波：a. 定义滤波器规格，计算其系数并绘制其传递函数。b. 使用设计好的滤波器对两个正弦波的混合信号进行滤波，并绘制输出信号。
a. 代码如下：
fs = 10e3;
fp = 1e3;
fc = 1.56e3;
ondp = .01;
ondc = .1;
ap = 1;
ac = 0;
[n, fo, mo, w] = remezord([fp fc], [ap ac], [ondp ondc], fs);
b = remez(n, fo, mo, w);
f1 = 1e3;
f2 = 1.56e3;
f1n = f1/fs;
f2n = f2/fs;
sig1 = sin(2*pi*f1n*[0:99]);
sig2 = sin(2*pi*f2n*[0:99]);
sig = sig1 + sig2;
y = filter(b, 1, sig);
[h, w] = freqz(b, 1);
fq = fs*w/(2*pi);
amp = 20*log10(abs(h));
figure;
plot(fq, amp);
axis([0 max(fq) min(amp) 10]);
xlabel('Frequency [Hz]');
ylabel('Amplitude [dB]');
grid;
b. 代码如下：
figure;
subplot(2, 2, 1);
plot(sig1);
title('First signal');
subplot(2, 2, 3);
plot(sig2);
title('Second signal');
subplot(2, 2, 2);
plot(sig);
title('Mixture of the two signals');
subplot(2, 2, 4);
plot(y);
title('Filtered signal');
30、二维有限脉冲响应（FIR）滤波器的设计方法与一维滤波器的设计方法类似。本练习的目的是说明图像的低通滤波。为此，考虑一个7×7像素的模板，中间的9个像素值为1，其他像素值为0。计算该滤波器的系数，并绘制相应的理想和实际传递函数。然后使用设计好的滤波器对MATLAB图像库中被方差为0.05的零均值白高斯噪声污染的“小丑”图像进行低通滤波。指出由于与一维滤波器相同的吉布斯现象导致的传递函数的波纹。
以下是MATLAB代码实现：
H = zeros(7,7);
H(3:5,3:5) = ones(3,3);
h = fsamp2(H);
figure;
subplot(2,2,1);
freqz2(h,32,32);
title('Actual transfer function');
[f1,f2] = freqspace([7 7]);
[x,y] = meshgrid(f1,f2);
subplot(2,2,2);
mesh(x,y,H);
title('Ideal transfer function');
load clown;
I = ind2gray(X,map);
I = imnoise(I,'gaussian',0,0.05);
subplot(2,2,3);
imshow(I);
title('Noisy image');
F = filter2(h,I);
subplot(2,2,4);
imshow(F);
title('Filtered image');
通过观察绘制的实际传递函数图，可以指出由于吉布斯现象导致的波纹。
31、一个包含20个二进制元素的数字信号在通信信道上传输。该信号具有以下特征： – 二进制速率：1000符号/秒， – 采样频率：10kHz， – 符号持续时间：1ms， – 编码：单极性不归零码。使用subplot函数绘制信号、匹配滤波器的脉冲响应及其输出信号。定义判决阈值和最佳判决时刻。
以下是MATLAB代码实现：
Rb = 1000; % 二进制速率
fs = 10000; % 采样频率
Ts = 1/fs; % 采样周期
Tb = 1/Rb; % 二进制元素持续时间
sequence=[0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1];
no = length(sequence); % 比特数
no_ech = no * Tb/Ts; % 样本数
t = [0:(no_ech-1)]*Ts; % 采样时刻
pulse = ones(1,fs/Rb); % 波形（单极性不归零）
x = (sequence'*pulse)'; % 信号生成
x = x(:);
t_time = Ts * [0:no_ech];
figure
subplot(311);
plot(t_time,[x;0]);
axis([0 0.021 0 1.2]);
grid
title('基带信号')
g = (pulse);
h = g(length(g):-1:1); % 匹配滤波器脉冲响应
subplot(312);
plot([0:size(h,2)+1]*Ts,[h 1 0]);
axis([0 0.021 0 1.2]);
grid
title('匹配滤波器脉冲响应')
y_fil = conv(h,x); % 匹配滤波器输出信号
y_fil = [y_fil(:)' 0];
subplot(313);
plot([0:size(y_fil,2)]*Ts,[0 y_fil]);
axis([0 0.021 0 1.2*max(y_fil)]);
grid
xlabel('时间 [s]');
title('匹配滤波器输出信号')
由于脉冲响应宽度为1ms，最佳判决时刻应为该量的倍数。对于等概率符号，根据最大似然准则，最佳判决阈值等于符号幅度的一半。匹配滤波器输出信号的幅度等于接收信号在匹配滤波器脉冲响应上投影的能量。
32、使用1000个符号的二进制序列和之前定义的信道生成眼图。注：这里假设之前定义的信道的相关参数（如脉冲形状、滤波器系数等）在代码中可正常使用。
以下为生成眼图的代码步骤：
sequence = round(rand(1, 1000));
x = (sequence' * pulse)'; % 信号生成
x = x(:);
y = filter(Bc, Ac, x); % 信号滤波
r = randn(size(y));
noise_power = 0.25;
y = y + sqrt(noise_power) * r; % 添加噪声
no_eyes = 2; % 眼图的眼数
no_block = fix(tot_ech / (no_eyes * no)); % 样本块的数量
index = [1:(no_eyes * no_ech)];
t_time = Ts * [0:no_eyes * no_ech];
time = Ts * [1:no_eyes * no_ech];
foo = zeros((no_eyes * no_ech), no_block - 1);
foo(:) = y_fil((no_eyes * no_ech) + 1: no_block * no_eyes * no_ech);
first_row = foo(no_eyes * no_ech, :);
first_row = [y_fil(maxindex) first_row(1:no_block - 2)];
foo = [first_row; foo];
plot(t_time, foo, 'b');
注：代码中
pulse
、
Bc
、
Ac
、
tot_ech
、
no
、
Ts
、
y_fil
、
maxindex
等变量需要根据实际情况进行定义。
33、本练习旨在比较三种二维滤波器（Prewitt、Sobel 和 Roberts）的结果。为此生成了一张包含垂直轮廓的图像，并为每种滤波器估计两个数值： – 检测概率：正确检测到的轮廓点的百分比； – 虚警概率：取自垂直轮廓外部但被视为轮廓点的点的百分比。为简化比较，设置检测阈值以获得大致相同的虚警概率水平。然后仅使用检测概率对这三种滤波器进行比较。
图像处理练习概述
练习目的
比较以下三种二维滤波器的结果：
Prewitt 滤波器
Sobel 滤波器
Roberts 滤波器
实验步骤
生成图像
：创建一个包含垂直轮廓的图像。
滤波器应用
：对每种滤波器进行应用。
性能评估
：
– 估计
检测概率
（正确检测轮廓点的百分比）
– 估计
虚警概率
（轮廓外被误判为轮廓点的百分比）
参数设置
设置
检测阈值
，使不同滤波器的虚警概率大致相同。
比较方法
在虚警概率一致的前提下，仅通过
检测概率
对三种滤波器进行比较。
34、考虑一个具有如下传递函数的系统，该传递函数假定未知，需要进行估计：1/(1 – 0.8z^(-1) – 0.5z^(-2))。编写一个 MATLAB 代码，使用递归最小二乘法（RLS）算法和包含 1023 个值的二进制伪随机序列（BPRS）来估计模型参数（此长度相对于系统的频率带宽而言足够大）。绘制估计系数随时间的变化曲线，并评估误差信号的自相关函数，该误差信号定义为系统实际系数对应的输出与估计系数对应的输出之间的差值。
以下是一个满足需求的 MATLAB 代码示例：
% 生成 BPRS 信号
n = 1023;
np = 3;
u = zeros(1, n);
u(1:10) = (-1).^(1:10);
for i = 11:n
u(i) = -u(i - 7) * u(i - 10); % 生成 BPRS
end
% 系统实际参数
para = [0.8, 0.5];
% 初始化变量
p = 1e8 * eye(np);
teta = zeros(np, 1);
v = randn(1, n);
y = zeros(1, n);
% 使用 RLS 算法估计参数
for i = 3:n
x = [y(i - 1), u(i - 1), u(i - 2)]';
y(i) = para * x + v(i);
k = (p * x) / (1 + x' * p * x);
teta = teta + k * (y(i) - x' * teta);
p = p - k * x' * p;
end
% 绘制估计系数随时间的变化曲线
figure;
for i = 1:np
plot(1:n, teta(i, :));
hold on;
end
axis([0 1000 -1.5 1.5]);
title('参数随时间的变化');
xlabel('时间');
grid on;
% 计算估计输出
ys = zeros(1, n);
for i = 3:n
ys(i) = [ys(i - 1), u(i - 1), u(i - 2)] * teta';
end
% 计算误差信号
err = ys - y;
average = mean(err);
err = err - average;
% 计算误差信号的自相关函数
cor_err = xcorr(err, err);
% 绘制误差信号的自相关函数
figure;
plot([-length(err) + 1:length(err) - 1], cor_err);
grid;
xlabel('时间');
title('误差信号的自相关函数');
此代码首先生成了一个长度为 1023 的 BPRS 信号，接着使用 RLS 算法来估计系统的参数。之后绘制了估计系数随时间的变化曲线，最后计算并绘制了误差信号的自相关函数。
© 版权声明
文章版权归作者所有，未经允许请勿转载。
相关文章
暂无评论...
