【生物医学信号处理及其MATLAB应用】Chapter2 数字信号处理基础

2.1 傅里叶变换及其意义

​ 成谐波关系的复指数函数在LTI系统中的分析：如果一个LTI系统的输入可以表示为周期复指数的线性组合，则输出也一定能表示成这种形式，并且输出线性组合中的加权系数与输入中对应的系数有关。

​ 傅里叶变换在LTI系统分析中的思想：将复杂的输入信号分解后分解成复指数信号的线性组合，那么系统的输出也能通过如上表所示的关系表达成相同复指数信号的线性组合，并且在输出中的每一个频率的复指数函数上乘以系统在那个频率的频率响应值。

​ 从表中时域和频域的关系还能得到如下规律：时域的离散必然导致频域的周期化；频域的离散必然导致时域的周期化。简单来说，就是一个域离散必然另外一个域是周期的；相反，如果一个域连续必然另外一个域是非周期的。

​ （本章仅对数字信号处理基础知识做简要介绍与回顾，以便为后续知识做铺垫，理论推导等具体内容请参考奥本海姆《信号与系统》和Sanjit《数字信号处理：基于计算机的方法》）

2.2 离散傅里叶变换

​ $DFT$变换对：

​ $$X(k)=DFT[x(n)]=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}kn},0\le k\le N-1$$（2-1）

​ $$x(n)=IDFT[X(k)]=\frac{1}{N}\sum _{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk},0\le n\le N-1$$（2-2）

​ 式（2-1）称为正变换，式（2-2）称为反变换。注意：这一对变换对中信号$x(n)$的长度为$N$，它的频谱$X(k)$点长也为$N$，则$x(n)$和$X(k)$具有唯一的映射对应关系。

​ 为了表示方便，一般用符号$W_N$来表示正交序列集中的基$e^{-j\frac{2\pi }{N}}$，即$W_N=e^{-j\frac{2\pi }{N}}$，因此，离散傅里叶变换对也可表示为：

​ $$X(k)=DFT[x(n)]=\sum _{n=0}^{N-1}x(n)W_N^{nk},0\le k\le N-1$$

​ $$x(n)=IDFT[X(k)]=\frac{1}{N}\sum _{k=0}^{N-1}X(k)W_N^{-nk},0\le n\le N-1$$

​ $W_N$具有下列性质：

​ 周期性：$W_N^n=W_N^{n+rN}$

​ 共轭对称性：$W_N^n=(W_N^{-n}{)}^{\ast }$

​ 可约性：$W_N^{rn}=W_{N/r}^n$或$W_{rN}^{rn}=W_N^n$

2.3 频谱分析和谱图表示

​ 1.幅度谱

​ 设$N$点长序列$x(n)$的$DFT$结果是$X(k)$，它是离散的复序列，可以用下式表示：

​ $$X(k)=|X(k)|e^{j\angle X(k)}$$

​ 带入傅里叶变换的反变换公式有：

​ $$x(n)=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{j[\frac{2\pi }{N}+\angle X(k)]}$$

​ 离散傅里叶变换的模$|X(k)|=\sqrt{[Re(X(k)){]}^2+[Im(X(k)){]}^2}$，表示信号$x(n)$的各复指数信号的频率分量的相对大小。

​ 2.相位谱

​ $\angle X(k)$表示相位角，$\angle X(k)=t{g}^{-1}(\frac{Im[X(k)]}{Re[X(k)]})$，它的大小不会影响各复指数频率分量的大小，但能提供这些频率的初始相位信息。$X(k)$对信号$x(n)$的性质有着显著的影响，因此一般包含了信号的大量信息，用相同的幅度谱和不同的相位谱得到的信号完全不同。

2.4 傅里叶变换的性质及其MATLAB验证

2.4.1 线性

x = rand(1,100);y = rand(1,100);
a = 2;b = 3;
X = fft(x);  % x的DFT
Y = fft(y);  % y的DFT
x_com = a*x + b*y;  % x和y的线性组合
X_com = fft(x_com);  % x_com的DFT

% 校验
X_check = a*X + b*Y;
error = max(abs(X_com - X_check));  % 差值
disp(error)

2.4.2 时间翻转特性：$$DFT[x(N-N)]=X(N-k)$$

x = 0:9;
N = length(x);
x1 = [x(1), flip(x(2:10))];  % 用flip函数把x翻转为0 9 8 7 6 5 4 3 2 1
X = fft(x);X1 = fft(x1);
X2 = [X(1),flip(X(2:10))];  % 用flip函数把X翻转
isequal(X1, X2)  % 比较X1与X2是否相等

2.4.3 序列的循环移位，使用circshift实现

x = 0:9;
f = circshift(x,8);  % circshift，使输出为f = 2 3 4 5 6 7 8 9 0 1

2.4.4 直观认识循环卷积和线性卷积的关系

% 线性卷积的实现
x = [1 1 1];y = [2 3 4 5];
e = conv(x, y);  % 求线性卷积
% 循环卷积的实现
x = [1 1 1];y = [2 3 4 5];N = 7;
X = fft(x, N);Y = fft(y, N);
f = ifft(X.*Y);  % 基于fft求循环卷积
f1 = cconv(x, y, N);  % cconv函数直接求循环卷积

2.4.5 共轭对称性

x = [0 1 2 3 4 4 3 2 1];
X = fft(x);
Xconj = fft(conj(x));
isequal(Xconj,X)

2.4.6 帕塞瓦尔（Parseval）定理

% 离散傅里叶变换的性质——帕塞瓦尔定理
% 信号的能量在时域和频域相等
% 离散形式 ∑u(k)² = 1/N( ∑U(k)U*(k) )
% Parseval's Relation
x = [0 1 2 3 4 5 6];
XF = fft(x);
% a代表着时域能量；b代表着频域能量
a = sum(x.*x)
b = round(sum(abs(XF).^2)/7)  % 用于舍入到最接近的整数

2.4.7 幅度谱和相位谱图及信号长度对谱图的影响

N = 100;
n = [0:1:(N-1)];  % 改变N的值，观察对幅度谱的影响，例如n取10个点，50个点，60个点，100个点，总结结果
x = cos(0.48*pi*n) + cos(0.52*pi*n);  % 画出x信号的幅度谱和相位谱
k = 1:1:N;
w = pi*k/N;
X = fft(x);
subplot(2,1,1);plot(w/pi, abs(X));
xlabel('以pi弧度为单位的频率');ylabel('幅度响应');
subplot(2,1,2);plot(w/pi, angle(X));
xlabel('以pi弧度为单位的相位');ylabel('相位响应');

2.5 频域分辨率

2.5.1 频域分辨率

​ 在时域中，数字信号的时间分辨率由采样间隔$\Delta t$决定，在频域中，频域分辨率$\Delta f$（两相邻谱线间的频率差值）由采样频率$f_s$和采样点数$N$决定：

​ $$\Delta f=\frac{f_s}{N}$$

​ 式中，$f_s$由采样定理决定，因此，要提高频域分辨率，就要增加采样点数$N$。如果数据量略有不足，传统方法是在数据尾部填0来解决，称为高密度频谱（The High Density Spectrum）。但是补零并不能提高频域分辨率，我们认为填入适当的现有数据会更好，称为高分辨率频谱（The High Resolution Spectrum）。

2.5.2 利用FFT做谱分析时的参数选择

​ 利用离散傅里叶变换进行频谱分析时会引起误差，常见的是由于离散傅里叶变换本身采样和截断过程所引起的混叠、泄漏、栅栏现象等。

​ 如果要分析的信号是连续信号，就必须先对信号采样，当采样频率比信号最高频率的两倍小，就会发生混叠现象，可以提高采样率来避免混叠现象；如果要分析的信号是周期信号，就必须对该信号截取一段来进行分析，即加了一个窗，便会发生泄漏现象，这时离散傅里叶变换的结果似乎是将原信号频谱扩展到整个频率范围，使得无法反映真正的频谱，要想减少泄漏可以通过加不同的窗函数来截取信号；离散傅里叶变换是对离散时间傅里叶变换的采样，它只给出频谱在离散点上的值，而无法反映这些点之间的频谱内容，这就是栅栏现象，改善栅栏现象的一种方法是信号后面补若干个零，通过补零来调整坐标上离散点的位置。

​ 通过上述讨论，利用FFT做谱分析时各参数的选择为：

​ ①采样频率应满足采样定理：$f_s\ge 2F_h$或$T\le 1/2F_h$

​ ②离散傅里叶变换的点长$N$为：$N=\frac{f_s}{\Delta f}\ge \frac{2F_h}{\Delta f}$，一般，为了FFT计算快速，$N$都尽量取成2的幂次

​ ③采集信号的持续时间为：$t_p=NT$

附录 信号多角度认知和滤波器的设计

①信号的时域和频域认识

%% 方波信号
clear, clc, close all
% 方波信号的产生
Fs = 1000;  % 采样率
t = 0:1/Fs:0.4;  % 样本点
N = length(t);  % 样本点个数
F = 30;  % 模拟频率为30Hz
y_sqr = square(2*pi*30*t);  % 方波信号的产生

% 方波信号时域绘制
figure;
plot(t, y_sqr);  % 绘制方波信号时域，设图中线宽1.5磅
title('方波信号时域图', 'fontsize', 20);
xlabel('时间/s', 'fontsize', 20);
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
ylim([-1.5, 1.5]);
grid minor

% 方波信号频域绘制
y_sqr_fft_abs = abs(fft(y_sqr));  % 幅度谱
y_sqr_fft_ang = angle(fft(y_sqr));  % 相位谱
f = (0:N-1)*Fs/N;  % Fs/N是频域采样率，f是频域的样本
figure;
plot(f, y_sqr_fft_abs, 'LineWidth', 1.5);  % 绘制幅度谱，设图中线宽1.5磅
title('方波信号幅度谱', 'fontsize', 20);
xlabel('频率/Hz', 'fontsize', 20);
ylabel('幅度', 'fontsize', 20);
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
grid minor
figure;
plot(f, y_sqr_fft_ang, 'LineWidth', 1.5);  % 绘制相位谱，设图中线宽1.5磅
title('方波信号相位谱', 'fontsize', 20);
xlabel('频率/Hz', 'fontsize', 20);
ylabel('幅度', 'fontsize', 20);
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
grid minor

%% chrip信号
clear, clc
Fs = 1000;
t = 0:1/Fs:1;
N = length(t);
y_chr = chirp(t, 0, 1, 150);  % chrip信号的产生

% chrip信号时域绘制
figure;plot(t, y_chr, 'linewidth', 1.5);
title('chrip信号时域图', 'fontsize', 20);
xlabel('时间/s', 'fontsize', 20)
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
ylim([-1.5,1.5]);
grid minor

% chrip信号频域绘制
y_chr_fft_abs = abs(fft(y_chr));  % 幅度谱
y_chr_fft_angle = angle(fft(y_chr));  % 相位谱
f = (0:N-1)*Fs/N;  % Fs/N是频域采样率，f是频域的样本
figure;plot(f, y_chr_fft_abs, 'LineWidth', 1.5);  % 幅度谱绘制
title('chrip信号幅度谱', 'fontsize', 20);
xlabel('频率/Hz', 'fontsize', 20);
ylabel('幅度', 'fontsize', 20);
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
grid minor
figure;
plot(f, y_chr_fft_angle, 'LineWidth', 1.5);  % 相位谱绘制
title('chrip信号相位谱', 'fontsize', 20);
xlabel('频率/Hz', 'fontsize', 20);
ylabel('相角/\pi', 'fontsize', 20);
set(gca, 'Fontsize', 20, 'LineWidth', 1.5);
grid minor

%% chrip信号分段
% 10段
numSeg = 10;  % 段数
len = floor(N/numSeg);
f_new = (0:len - 1)*Fs/len;
figure;
for i = 1:numSeg
    y_seg = y_chr((i - 1)*len + 1:i*len);
    y_seg_fft = abs(fft(y_seg));
    subplot(3,4,i);
    plot(f_new, y_seg_fft, 'LineWidth', 1.5);
end

% 50段
numSeg2 = 50;
len = floor(N/numSeg2);
f_new = (0:len - 1)*Fs/len;
figure;
for i = 1:numSeg2
    y_seg = y_chr((i - 1)*len + 1:i*len);
    y_seg_fft = abs(fft(y_seg));
    subplot(6,9,i);
    plot(f_new, y_seg_fft, 'LineWidth', 1.5);
end

% 100段
numSeg3 = 100;
len = floor(N/numSeg3);
f_new = (0:len - 1)*Fs/len;
figure;
for i = 1:numSeg3
    y_seg = y_chr((i - 1)*len + 1:i*len);
    y_seg_fft = abs(fft(y_seg));
    subplot(10,10,i);
    plot(f_new, y_seg_fft, 'LineWidth', 1.5);
end

②滤波器的相位谱特性对滤波结果的影响

% 构造一个任意离散信号x
fs = 100;
t = 0:1/fs:10;
x = 1.5*sin(2*pi*t) + 2*cos(3*pi*t) + exp(0.3*t);
n = length(x);
% 计算DFT频谱
x_fft = fft(x);
% 改造相位信息
x_fft_abs = abs(x_fft);
x_fft_angle = angle(x_fft);
x1_angle = x_fft_angle - exp(5*t);  % 任意相位
x2_angle = x_fft_angle + 8*t;  % 线性相位
% 对改造后的频谱做IDFT
x1 = ifft(x_fft_abs.*exp(1i*x1_angle), 'symmetric');
x2 = ifft(x_fft_abs.*exp(1i*x2_angle), 'symmetric');
figure;plot(t, x, t, x1, t, x2);
legend('原始信号', '叠加任意相位后的信号', '叠加线性相位后的信号');

③IIR数字高通滤波器的设计

clear;
fs = 2000;
t = 0:1/fs:1;
y1 = 0 + sqrt(1)*randn(1, length(t));  % 滤波前信号
IIR = Butterworth;
y2 = filter(IIR, y1);
n1 = length(abs(fft(y1)));
n2 = length(abs(fft(y2)));
f1 = (0:n1 - 1)*fs/n1;
f2 = (0:n2 - 1)*fs/n2;
subplot(3, 1, 1);
plot(t, y1, t, y2);title('滤波前后信号对比');legend('滤波前信号', '滤波后信号');
subplot(3, 1, 2);
plot(f1, abs(fft(y1)), f2, abs(fft(y2)));title('幅度谱');legend('原始信号', '滤波后信号');
axis([1, 1000, 0, 100]);
xlabel('频率/Hz');ylabel('幅度');
subplot(3, 1, 3);
plot(angle(fft(y1)));hold on;
plot(angle(fft(y2)));
title('相位谱');legend('原始信号', '滤波后信号');
axis([0, 1000, -5, 5]);
ylabel('相位');

④FIR滤波器设计

clc;
fs = 1000;
t = 0:1/fs:1;
y1 = 5 + sqrt(0.2)*randn(1, length(t));  % 滤波前信号
FIR = Hamming;
y2 = filter(FIR, y1);
n1 = length(abs(fft(y1)));
n2 = length(abs(fft(y2)));
f1 = (0:n1 - 1)*fs/n1;
f2 = (0:n2 - 1)*fs/n2;
subplot(3, 1, 1);
plot(t, y1, t, y2);title('滤波前后信号对比');legend('滤波前信号', '滤波后信号');
subplot(3, 1, 2);
plot(f1, abs(fft(y1)), f2, abs(fft(y2)));title('幅度谱');legend('原始信号', '滤波后信号');
axis([1, 1000, 0, 100]);
xlabel('频率/Hz');ylabel('幅度');
subplot(3, 1, 3);
plot(angle(fft(y1)));hold on;
plot(angle(fft(y2)));
title('相位谱');legend('原始信号', '滤波后信号');
axis([0, 1000, -5, 5]);
ylabel('相位');

⑤滤波和卷积的关系

clear all;
wp = 0.67*pi;
ws = 0.53*pi;
M = 3.32*pi/(wp - ws);
N = ceil(2*M + 1);
hamwin = hamming(N + 1);  % 加汉明窗
wc = (wp + ws)/2/pi;
b = fir1(N, wc, hamwin);

Ft = 2000;
x = randn(1, 0.2*Ft + 1);
x_filter = filter(b, 1, x);  % 滤波
x_conv = conv(b, x);  % 卷积
plot(1:401, x_filter, 'o-');hold on
plot(1:450, x_conv, '*-');
set(gcf, 'color', 'white');
legend('滤波后的信号', '卷积后的信号');
title('滤波和卷积的比较');