MATLAB编程与信号处理实践

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;
}

1、定义一个4×3的矩阵，除第一行元素为1外，其余元素均为0。

b = ones (1,3); m = zeros (4,3); m(1,:) = b

2、考虑向量对 (x1, y1) 和 (x2, y2)。定义向量 x，使得：当 y1(j) < y2(j) 时，x(j) = 0；当 y1(j) = y2(j) 时，x(j) = x1(j)；当 y1(j) > y2(j) 时，x(j) = x2(j)。请编写 MATLAB 代码实现该功能。

可以使用 MATLAB 代码来实现该功能，以下是实现代码：

% 假设 x1, y1, x2, y2 已经定义
x = zeros(size(y1)); % 初始化 x 为全零向量，长度与 y1 相同

idx_equal = y1 == y2;  % 找到 y1 等于 y2 的索引
idx_greater = y1 > y2; % 找到 y1 大于 y2 的索引

x(idx_equal) = x1(idx_equal);   % 在 y1 等于 y2 的位置，将 x1 对应位置的值赋给 x
x(idx_greater) = x2(idx_greater); % 在 y1 大于 y2 的位置，将 x2 对应位置的值赋给 x

3、生成并绘制信号：对于 0 ≤ t ≤ 2，以 0.01 为增量，y(t) = sin(2πt)，然后使用 decimate 函数分别以 2 和 16 为下采样因子对其进行下采样。

t = 0:0.01:2;
y = sin(2*pi*t);
subplot(311);
plot(t,y);
ylabel('sin(2*pi*t)');
title('Original signal');

t2 = decimate(t,2);
t16 = decimate(t2,8);
y2 = decimate(y,2);
y16 = decimate(y2,8);

subplot(312);
plot(t2,y2);
ylabel('sin(2*pi*t)');
title('Undersampled signal with a factor 2');

subplot(313);
plot(t16,y16);
ylabel('sin(2*pi*t)');
xlabel('Time t');
title('Undersampled signal with a factor 16');

print -eps file_name

4、在极坐标中绘制具有以下传递函数的滤波器的极点：$H(z)=frac{1}{1 + az^{-1}+bz^{-2}}$，其中$a$和$b$的值由用户输入，函数返回极点。（使用

roots

和

polar

命令）

以下是实现该功能的 MATLAB 代码：

function c = filter_bis(a,b)
    c = roots([1 a b]);
    fprintf('The poles are:
');
    z1 = c(1,:);
    z2 = c(2,:);
    if (abs(z1) > 1 | abs(z2) > 1)
        fprintf('There is at least an instable pole
');
    else
        clf;
        figure;
        polar(angle(z1), abs(z1), 'r+');
        hold on;
        polar(angle(z2), abs(z2), 'r+');
        legend('Polar plot of H(z) poles', 0);
    end
end

5、生成信号：x(t) = A⋅sin(2πft + φ) + b(t)，其中 t 取值范围是 0 到 1024，φ 是在[0, 2π]上均匀分布的随机变量，b(t) 是均值为零、方差为 1 的高斯白噪声（使用 rand 和 randn 函数）。A 和 f 由用户选择。使用周期图法和相关图法（使用 fft 和 fftshift 函数）估计 x(t) 的均值、自相关函数（xcorr）和频谱。将得到的结果与理论结果进行比较。通过改变 A 来控制信噪比。

以下是实现该功能的 MATLAB 代码：

function noisy_sin(A,f)
    N = 1024; % 计算频率的数量
    nech = 1024; % 样本数量
    t = 0:nech;
    b = randn(1,nech+1);
    phi = 2*pi*rand(1);
    x = A*sin(2*pi*f*t+phi)+b;
    fprintf('Mean of x(t) = %2.4f',mean(x));
    fprintf('Mean of b(t) = %2.4f',mean(b));
    cx = xcorr(x);
    % 基于相关图的频谱估计:
    sx_correlo = (abs(fft(cx,N))).^2;
    sx_correlo = sx_correlo / max(sx_correlo);
    % 前 N/2 个值对应 0<=f<0.5
    % 后 N/2 个值对应 0.5<=f<1 (或 -0.5<=f< 1)
    sx_correlo = fftshift(sx_correlo);
    % 频谱以 0 为中心:
    % 前 N/2 个值对应 -0.5<=f<0
    % 后 N/2 个值对应 0<=f<0.5
    % 基于周期图的频谱估计:
    sx_periodo = (abs(fft(x,N))).^2;
    sx_periodo = sx_periodo / max(sx_periodo);
    sx_periodo = fftshift(sx_periodo);
    % 带噪正弦信号的 SNR 估计
end

要比较结果与理论结果，需要根据理论公式手动计算均值、自相关函数和频谱，然后与代码输出的结果进行比较。改变 A 的值可以控制信噪比，因为 A 是正弦信号的幅度，A 越大，信号越强，信噪比越高。

6、编写一个MATLAB代码来验证帕塞瓦尔定理：对于信号x[n]，当n取值为0到7时，x[n]=n；n取其他值时，x[n]=0。验证公式sum_{n = 0}^{N – 1}|x[n]|^2 = (1/N) * sum_{k = 0}^{N – 1}|X[k]|^2 。

以下是用于验证该定理的MATLAB代码：

% 定义信号x[n]
x = 0:7;
N = length(x);

% 计算信号的平方和
sum_x_squared = sum(abs(x).^2);

% 计算信号的离散傅里叶变换
X = fft(x);

% 计算DFT系数的平方和
sum_X_squared = sum(abs(X).^2);

% 验证帕塞瓦尔定理
result = sum_x_squared - (1/N)*sum_X_squared;
if abs(result) < 1e-10
    disp('帕塞瓦尔定理验证通过');
else
    disp('帕塞瓦尔定理验证失败');
end

在上述代码中，首先定义了信号

x[n]

，然后计算信号的平方和

sum_x_squared

。接着，使用

fft

函数计算信号的离散傅里叶变换

X

，并计算DFT系数的平方和

sum_X_squared

。最后，通过比较等式两边的值来验证帕塞瓦尔定理。如果差值小于一个很小的数（如

1e-10

），则认为定理验证通过。

7、考虑两个独立同分布（i.i.d.）的随机变量X和Y，其概率密度函数（pdf）是参数λ = 2的指数函数。确定Z = X + Y的分布，将得到的分布与解析分布进行比较。然后计算对应于Z和 -X的概率密度函数的卷积。已知Z + (-X) = Y，如何解释其与Y分布的差异？

Z的分布通过对X和Y的概率密度函数进行卷积得到，即

Pz = conv(Px, Px) * dx

；

解析分布为

Pz_analytic = lambda^2 * exp(-lambda * x) .* x

。

对应于Z和 -X的概率密度函数的卷积为

Czx = conv(Pz, Px(length(Px):-1:1)) * dx

。

虽然

Z + (-X) = Y

，但由于 -X和Z不独立（

Z = X + Y

），所以Z和 -X的概率密度函数的卷积不等于Y的概率密度函数，不能通过两个概率密度函数的卷积来得到和的概率密度函数。

8、1000名年龄在20岁至50岁之间的人参加了一项统计调查。他们被要求提供自己的身高和体重，并将这些数据与标准值进行比较。实际值与标准值之间的差异是一个随机高斯变量，其方差取决于所考虑的参数。估算每对变量（体重，身高）、（年龄，体重）和（年龄，身高）的相关系数。

（体重，身高）的相关系数为0.8466；

（年龄，体重）的相关系数为0.2655；

（年龄，身高）的相关系数为0.0155。

9、考虑一个均值为 1、方差为 0.25 的高斯白随机过程 P 的 2048 个样本。使用不同的测试方法来检验它是否为高斯随机变量。

随机过程生成：

matlab
P = 0.5 * randn(1, 2048) + 1;


– 计算均值：

matlab
mp = mean(P);


– 计算方差：

matlab
vp = var(P);


第一种测试：

使用最简单的测试方法，即随机过程直方图。

– 设置

bins = 20

，并绘制直方图：

matlab
hist(P, bins);


– 从得到的直方图可以很容易看出，它与高斯概率密度函数非常相似。

10、一家工厂供应一系列每批1000件的特定类型产品。随机选取了一些批次进行检查，有k个次品的被检查批次数量nk如下：当k为0时，nk为7；当k为1时，nk为12；当k为2时，nk为23；当k为3时，nk为27；当k为4时，nk为41；当k为5时，nk为36；当k为6时，nk为25；当k为7时，nk为20；当k为8时，nk为9；当k为9时，nk为5；当k为10时，nk为2；当k为11时，nk为1；当k为12时，nk为2。使用卡方（Chi2）检验来验证或拒绝每批次次品数量服从泊松分布的假设。

通过计算卡方值并与临界值比较来判断每批次次品数量是否服从泊松分布。若卡方值小于临界值，则认为随机变量服从泊松分布；若卡方值大于等于临界值，则认为随机变量不服从泊松分布。具体判断结果需运行代码得出。

11、源分离是信号处理领域的一个新问题，它是利用q个混合信号来恢复p个相互独立的未知源，这p个源用向量S(t) = {si(t)}表示，其中1 ≤ i ≤ p，设Y(t)是观测向量，其分量为yi(t)，1 ≤ i ≤ q，它们是p个源的未知函数。本练习强调高阶统计量在进行盲源或半盲源分离中的重要性，还需考虑一个未知源向量S和一个观测向量Y（源信号的混合）。请简述源分离问题包含的主要内容。

源分离问题包含的主要内容有：

定义

：利用混合信号恢复独立未知源

相关向量表示

：

源向量 $ S(t) $

观测向量 $ Y(t) $

练习目的

：强调高阶统计量对盲源或半盲源分离的重要性

需考虑的要素

：

未知源向量 $ S $

观测向量 $ Y $

12、计算并绘制一维数字信号的频谱，该信号定义为：当n在0到7之间时，x[n] = 1；其他情况，x[n] = 0。

可以使用以下 MATLAB 代码进行计算和绘图：

x = ones(1,8);
X = fft(x,128);
subplot(3,1,1);
stem(x);
xlabel('n');
ylabel('x[n]');
subplot(3,1,2);
stem(abs(X));
xlabel('k');
ylabel('abs(X[k])');
subplot(3,1,3);
stem(angle(X));
xlabel('k');
ylabel('angle(X)');

13、编写一个MATLAB代码，计算MATLAB中“Clown”和“Gatlin2”两幅图像的二维傅里叶变换。绘制相应的频率表示，并通过逆变换重建这两幅图像，但交换它们的相位谱。

以下是实现该功能的MATLAB代码：

clear; clf
%%%% IMAGE LOADING
% Loading the image x from the file clown.mat
load clown; x = X;
% Loading the image y from the file gatlin2.mat
load gatlin2; y = X;
% Resizing x and y to have the same size
l = min(size(x,1),size(y,1));
c = min(size(x,2),size(y,2));
x = x(1:l,1:c);
y = y(1:l,1:c);
%%%%% RECONSTRUCTION OF THE NEW IMAGES
% Calculation of the FT of x and y
X = fft2(x);
Y = fft2(y);
% Reconstruction of z1 using the magnitude of X and the phase of Y
z1 = real(ifft2(abs(X).*exp(1i*angle(Y))));
% Reconstruction of z2 using the magnitude of Y and the phase of X
z2 = real(ifft2(abs(Y).*exp(1i*angle(X))));
%%%%%% DISPLAYING THE IMAGES
% Displaying the FT of x and y: magnitude and phase
fh = ([1:l]-l/2)/l;
fv = ([1:c]-c/2)/c;
figure(1);
subplot(3,2,1);
image(x);
axis('image');
title('image x');
axis off
subplot(3,2,2);
image(y);
axis('image');
title('image y');
axis off
subplot(3,2,3);
imagesc(fh,fv,log(abs(fftshift(X.^2))));
axis('image');
title('Magnitude of X');
xlabel('horizontal frequencies');
ylabel('vertical frequencies')
subplot(3,2,4);
imagesc(fh,fv,log(abs(fftshift(Y.^2))));
axis('image');
title('Magnitude of Y');
xlabel('horizontal frequencies');
ylabel('vertical frequencies')
subplot(3,2,5);
imagesc(fh,fv,angle(fftshift(X)));
axis('image');
title('X phase');
xlabel('horizontal frequencies');
ylabel('vertical frequencies')
subplot(3,2,6);
imagesc(fh,fv,angle(fftshift(Y)));
axis('image');
title('Y phase');

14、编写一个MATLAB代码，用于说明信号的观测时间对其频谱的影响。考虑一个10 Hz的正弦信号，采样频率为100 Hz，并处理以下情况： – 信号周期数完整； – 在计算频谱之前对信号进行零填充； – 信号周期数不完整。

以下是处理不同情况的MATLAB代码示例：

情况1：信号周期数完整

% 时间支持定义（50个样本 = 5个周期）
t = 0:0.01:0.5-0.01;
x = sin(2*pi*10*t); % 10 Hz正弦信号
N = length(x);
X = fft(x,N);
fp = (0:N-1)/N/0.01;
fp = fp - 1/2/0.01;
stem(fp,abs(fftshift(X)));
axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);
xlabel('Frequency [Hz]')
ylabel('Amplitude spectrum')

情况2：信号零填充

% 时间支持定义（50个样本 = 5个周期）
t = 0:0.01:0.5-0.01;
x = sin(2*pi*10*t); % 10 Hz正弦信号
N = length(x);
N_padded = 2*N; % 零填充到原来长度的2倍
X_padded = fft([x zeros(1,N_padded - N)],N_padded);
fp_padded = (0:N_padded-1)/N_padded/0.01;
fp_padded = fp_padded - 1/2/0.01;
stem(fp_padded,abs(fftshift(X_padded)));
axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);
xlabel('Frequency [Hz]')
ylabel('Amplitude spectrum')

情况3：信号周期数不完整

% 时间支持定义（45个样本 = 4.5个信号周期）
t = 0:0.01:0.5-0.01-1/20;
x = sin(2*pi*10*t); % 10 Hz正弦信号
N = length(x);
X = fft(x,N);
fp = (0:N-1)/N/0.01;
fp = fp - 1/2/0.01;
stem(fp,abs(fftshift(X)));
axis([ -1/(2*0.01) 1/(2*0.01) 0 25]);
xlabel('Frequency [Hz]')
ylabel('Amplitude spectrum')
axis([-50 50 0 20])

15、考虑一个由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')

16、计算数字信号 $x[n]=n + 50cos(frac{2pi 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

。

17、以下矩阵：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(log(1 + abs(Xfsh))); 
title('DFT2D of the image');

subplot(223); 
imagesc(log(1 + abs(Xc))); 
title('DCT2D of the image');

18、考虑一个幅值为 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。

t = (1:128);
f = 12.5; % 信号频率

% 情况 a: Fs = 64 Hz，Nfft = 128
Fs1 = 64;
y1 = sin(2*pi*f/Fs1*t);
sig1 = y1;
Nfft1 = 128;
y1_fft = abs(fft(sig1, Nfft1));
f_axis1 = [0:Fs1/Nfft1:(Fs1 - Fs1/Nfft1)];

% 情况 b: Fs = 128 Hz，Nfft = 128
Fs2 = 128;
y2 = sin(2*pi*f/Fs2*t);
sig2 = y2;
Nfft2 = 128;
y2_fft = abs(fftshift(fft(sig2, Nfft2)));
f_axis2 = [-Fs2/2:Fs2/Nfft2:(Fs2/2 - Fs2/Nfft2)];

% 情况 c: Fs = 128 Hz，Nfft = 256
Nfft3 = 256;
y3_fft = abs(fftshift(fft(sig2, Nfft3)));
f_axis3 = [-Fs2/2:Fs2/Nfft3:(Fs2/2 - Fs2/Nfft3)];

% 绘制频谱图
subplot(3,1,1);
stem(f_axis1, y1_fft);
grid on;
axis([0 64 0 80]);
title('Fs = 64 Hz, Nfft = 128 points');
ylabel('Spectral amplitude');

subplot(3,1,2);
stem(f_axis2(Nfft2/2 + 1:Nfft2), y2_fft(Nfft2/2 + 1:Nfft2));
grid on;
title('Fs = 128 Hz, Nfft = 128 points');
ylabel('Spectral amplitude');

subplot(3,1,3);
stem(f_axis3, y3_fft);
grid on;
title('Fs = 128 Hz, Nfft = 256 points');
xlabel('Frequency (Hz)');
ylabel('Spectral amplitude');