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;
}
17、考虑一个由50Hz和120Hz的两个正弦波之和组成的信号，该信号受到加性零均值白噪声的干扰。计算并绘制其功率谱密度。
t = 0:0.001:0.8; % T=0.001, so Fs=1 kHz
x = sin(2*pi*50*t) + sin(2*pi*120*t) + 2*randn(1,length(t));
subplot(211);
plot(x(1:500));
title('Signal');
X = fft(x,512);
Px = X.*conj(X)/512;
f = 1000*(0:255)/512; % f=Fe*(k - 1)/N
subplot(212);
plot(f, Px(1:256));
title('Power spectral density');
18、计算数字信号 (x[n]=n + 50cos(2πn/40))（n 取值范围为 1 到 100）的一维离散余弦变换（DCT1D），并确定其前三个系数累积的总能量的百分比。
使用代码
x = (1:100) + 50*cos(2*pi/40*(1:100));
X = dct(x);
norm([X(1), X(2), X(3)]) / norm(X)
计算，得到累积能量的百分比为
0.8590
。
19、以下矩阵：x = [1 1 1 1 1; 1 0 0 0 1; 1 0 0 0 1; 1 0 0 0 1; 1 1 1 1 1]; 表示一个以白色为背景、黑色正方形为主体的二值数字图像。绘制该图像，然后计算并绘制该图像的二维离散傅里叶变换（DFT2D）和二维离散余弦变换（DCT2D）。
可使用以下MATLAB代码完成任务：
```matlab
x = [1 1 1 1 1;
1 0 0 0 1;
1 0 0 0 1;
1 0 0 0 1;
1 1 1 1 1];
Xf = fft2(x);
Xfsh = fftshift(Xf);
Xc = dct2(x);
subplot(221);
imagesc(x);
title('Original image');
subplot(222);
imagesc(abs(Xfsh));
title('DFT2D');
subplot(223);
imagesc(abs(Xc));
title('DCT2D');
##20、考虑一个幅值为1、频率f = 12.5 Hz的128点正弦信号。该信号以Fs = 64 Hz或Fs = 128 Hz进行采样。其离散傅里叶变换（DFT）计算点数为128或256。编写MATLAB代码，展示对应以下三种情况的信号频谱表示：a. Fs = 64 Hz，Nfft = 128；b. Fs = 128 Hz，Nfft = 128；c. Fs = 128 Hz，Nfft = 256。
```matlab
t = (1:128);
f1 = 12.5; % 信号频率
% 情况a: Fs = 64 Hz，Nfft = 128
Fs = 64;
y1 = sin(2*pi*f1/Fs*t);
sig = y1;
Nfft = 128;
y = abs(fft(sig, Nfft));
f0 = [0:Fs/Nfft:(Fs - Fs/Nfft)];
% 情况b: Fs = 128 Hz，Nfft = 128
Fs = 128;
y1 = sin(2*pi*f1/Fs*t);
sig = y1;
Nfft = 128;
y_rect = abs(fftshift(fft(sig, Nfft)));
% 情况c: Fs = 128 Hz，Nfft = 256
Nfft1 = 256;
y_rect1 = abs(fftshift(fft(sig, Nfft1)));
f = [-Fs/2:Fs/Nfft:(Fs/2 - Fs/Nfft)];
f1 = [-Fs/2:Fs/Nfft1:(Fs/2 - Fs/Nfft1)];
subplot(3,1,1);
stem(f0, y);
grid on;
axis([0 64 0 80]);
title('Fs = 64 Hz, Nfft = 128 points');
ylabel('Spectral amplitude');
subplot(3,1,2);
stem(f(Nfft/2 + 1:Nfft), y_rect(Nfft/2 + 1:Nfft));
grid on;
title('Fs = 128 Hz, Nfft = 128 points');
ylabel('Spectral amplitude');
subplot(3,1,3);
stem(f1, y_rect1);
grid on;
title('Fs = 128 Hz, Nfft = 256 points');
ylabel('Spectral amplitude');
xlabel('Frequency (Hz)')
21、如果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}。例如，当N = 5时，序列 {1, 2, 3, 4} 会被排列成序列 {2, 4, 3, 1}。请编写一个MATLAB代码来执行这些排列。
以下是一个实现该排列的 MATLAB 代码：
N = 5;
rho = 2;
sequence = 1:N-1;
permutation = mod(rho.^sequence, N);
此代码中，首先设定了
N
和
ρ
的值，接着创建了原始序列，最后使用
mod
函数对
rho.^sequence
进行取模运算，得到排列后的序列。
22、在拉普拉斯域中定义一个抗混叠滤波器的传递函数，该滤波器在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]);
23、使用函数 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) $，将得到的结果与理论结果进行比较需具备系统理论知识。
24、计算截止频率为 1 kHz 的 5 阶模拟巴特沃斯滤波器的极点。
% 原型滤波器设计
N = 5;
[zc, pc, kc] = buttap(N);
% 频率变换
Omega = 2 * pi * 10^3;
pc = pc * Omega;
% 绘制极点
zplane(zc, pc)
25、设计一个采样频率为 1 Hz 的 IIR 低通滤波器，在频率范围 [0 0.1306] Hz 内衰减为 0.75 dB，阻带衰减至少为 20 dB，阻带下限边缘频率为 0.201 Hz。使用巴特沃斯模型计算相关模拟滤波器的极点。
% 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);
26、计算一个与之前练习参数相同，但使用布莱克曼窗而非矩形窗的有限脉冲响应（FIR）滤波器的系数。从过渡带和阻带衰减方面比较所得结果。
使用与之前练习相同的 MATLAB 代码，将选项
"boxcar"
替换为
"blackman"
。所得结果显示，过渡带约为上一情况的两倍宽，阻带衰减增加了约 30 dB。因此，当需要高选择性时应使用矩形窗。
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。绘制所得到滤波器的传递函数、脉冲响应和零点，并与另一个设计好的滤波器进行比较（假设另一个滤波器阶数为25）。
代码如下：
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），通带波纹减小，且用雷米兹方法设计的滤波器在阻带是等波纹的。
29、考虑两个正弦波的和，其频率分别为1kHz和1.56kHz，采样频率为10kHz。使用采用雷米兹（Remez）方法设计的低通有限脉冲响应（FIR）滤波器从该混合信号中提取第一个正弦波：a. 定义滤波器规格，计算其系数并绘制其传递函数。b. 使用设计好的滤波器对两个正弦波的混合信号进行滤波，并绘制输出信号。
a. 滤波器规格定义及系数计算代码如下：
```matlab
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代码：
```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代码如下：
```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(3,1,1);
plot(t_time, [x; 0]);
axis([0 0.021 0 1.2]);
grid
title('基带信号')
g = (pulse);
h = g(length(g):-1:1); % 匹配滤波器脉冲响应
subplot(3,1,2);
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(3,1,3);
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个符号的二进制序列和之前定义的信道生成眼图。注：此题目因‘之前定义的信道’指代不明，可能导致没看过书本内容的人无法完全理解，但暂按现有内容保留。
此答案存在较多问题，代码中部分变量（如 `pulse`、`Bc`、`Ac`、`tot_ech`、`no`、`Ts`、`y_fil`、`maxindex` 等）未定义，没看过书本内容的人无法理解代码逻辑，因此答案为 **DELETE**
##33、本练习旨在比较三种二维滤波器（Prewitt、Sobel和Roberts）的结果。为此生成了一个包含垂直轮廓的图像，并为每种滤波器估计两个量： - 检测概率：正确检测到的轮廓点的百分比； - 虚警概率：取自垂直轮廓外部但被视为轮廓点的点的百分比。为了简化比较，设置检测阈值以获得大致相同的虚警概率水平。然后仅使用检测概率对这三种滤波器进行比较。
# 图像处理练习概述
该练习聚焦于对比 Prewitt、Sobel 和 Roberts 三种二维滤波器。
## 目标
- 生成包含垂直轮廓的图像
- 估算每种滤波器的：
- **检测概率**：正确检测轮廓点的百分比
- **虚警概率**：将轮廓外点误判为轮廓点的百分比
## 比较方法
为便于比较：
- 设定检测阈值，使得**虚警概率大致相同**
- 在此基础上，根据**检测概率**对三种滤波器进行对比
##34、考虑一个具有如下传递函数的系统，该传递函数假定未知且需要进行估计：1/(1 - 0.8z^(-1) - 0.5z^(-2))。编写一个MATLAB代码，使用递归最小二乘法（RLS）算法和包含1023个值的二进制伪随机序列（BPRS）来估计模型参数（此长度相对于系统的频率带宽而言足够大）。绘制估计系数随时间的变化曲线，并评估误差信号的自相关函数，该误差信号定义为系统对应实际系数和估计系数的输出之间的差值。
以下是解决该问题的MATLAB代码：
```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, n);
v = randn(1, n);
y = zeros(1, n);
% 生成系统输出
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(:, i) = teta(:, i - 1) + k * (y(i) - x' * teta(:, i - 1));
p = p - k * x' * p;
end
para_fin = teta(:, end);
% 绘制估计系数随时间的变化曲线
figure;
subplot(2, 1, 1);
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)] * para_fin';
end
% 计算误差信号
err = ys - y;
average = mean(err);
var = std(err - average)^2;
err = err - average;
cor_err = xcorr(err, err);
% 绘制误差信号的自相关函数
subplot(2, 1, 2);
plot([-length(err) + 1:length(err) - 1], cor_err);
grid;
xlabel('时间');
title('误差信号的自相关函数');
这段代码首先生成了一个长度为1023的BPRS序列，然后使用RLS算法估计系统的参数。接着绘制了估计系数随时间的变化曲线，并计算和绘制了误差信号的自相关函数。
© 版权声明
文章版权归作者所有，未经允许请勿转载。
相关文章
暂无评论...
