《数字信号处理》大作业

《数字信号处理》大作业

目录

一、题目分析... 2

1.1 HWDATA.MAT文件... 2

1.2 程序设计思路... 2

1.3 FFT算法原理... 2

1.4 Butterworth滤波器设计... 3

二、系统框图：... 4

三、代码实现和运行结果... 4

3.1 导入数据... 4

3.2绘制原始信号时域波形... 4

3.3绘制原始信号的频域波形... 5

3.4设计Butterworth带通滤波器... 6

3.5绘制目标信号的时域波形... 6

3.5绘制目标信号的频域波形... 7

四、总结... 7

程序源码：

创建.m文件，将下述代码复制到.m文件当中，即可在matlab中运行

load('HWDATA.mat');
%导入HWDATA.mat后，产生的变量名叫 s1，为了方便使用，这里把s1重命名为HW
HW=s1;
clear s1;

fs = 500; % 采样频率为500Hz
t = (0:length(HW)-1)/fs; % 生成时间序列

figure;
subplot(4,1,1);
plot(t, real(HW));
ylabel('幅值');
xlabel('时间(s)');
title('原始信号');

N = length(HW); % 采样点数
Y = fft(HW); % 进行傅里叶变换
f = (0:N-1)*(fs/N); % 计算出频率坐标轴
P = abs(Y/N); % 计算出功率谱

%figure;
subplot(4,1,2);
plot(f, P);
xlim([0, fs/2]);
ylabel('功率');
xlabel('频率(Hz)');
%title('原始信号频谱');

% 设计Butterworth滤波器

Low = 94;   % 截止频率下限
High = 96;  % 截止频率上限
order = 4;   % 阶数为4

% 计算滤波器系数
[b,a] = butter(order, [2*Low/fs, 2*High/fs], 'bandpass');  

% 对信号进行滤波
HW_filtered = filtfilt(b, a, HW);

% 绘制滤波后的信号和频谱
subplot(4,1,3);
plot(t,real(HW_filtered));
ylabel('幅值');
xlabel('时间(s)');
title('滤波后的信号');

subplot(4,1,4);
NFFT = 2^nextpow2(length(HW_filtered));
Y = fft(HW_filtered, NFFT)/length(HW_filtered);
f = fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(Y(1:NFFT/2+1)));
xlim([90 100]);
ylabel('功率');
xlabel('频率(Hz)');

一、题目分析

已知某组数字信号（见大作业数据压缩包中 HWDATA.mat 文件），该信号中除了目标信号之外还掺杂有强噪声，但噪声与目标信号的频率不重叠，要求采用本学期已学的知识对该信号进行处理，并得到无噪声的目标信号的频谱图。（20 分）

1.1HWDATA.MAT文件

HWDATA.MAT文件中的数据是由500个采样点构成的时域信号。

1.2程序设计思路

因为目标信号与噪音信号的频率不重叠，所以可以通过观察原始信号的频域波形，来确定目标信号的频率，再设计符合条件的滤波器处理原始信号，就可以的到目标信号了。

1.3FFT算法原理

信号由时域到频域的变化可以使用FFT算法。FFT是一种DFT的高效算法，称为快速傅立叶变换（fast Fourier transform）。傅里叶变换是时域到频域变换分析中最基本的方法之一。FFT算法可分为按时间抽取算法和按频率抽取算法。

   DFT对于X（K）的每个K值，需要进行4N次实数相乘和（4N-2）次相加，对于N个k值，共需N x N乘和N（4N-2）次实数相加。FFT改变了运算路径，减小它的运算量。FFT利用DFT中的周期性和对称性，使整个DFT的计算变成一系列迭代运算，可大幅度提高运算过程和运算量。

   FFT分为时间抽取和频率抽取，最常用的时域抽选方法是按奇偶将长序列不断地变为短序列，结果使输入序列为倒序，输出序列为顺序排列。

   对于时域抽取，设N点序列x(n),将x(n)按奇偶分组，一个N点DFT分解为两个 N/2点的DFT，由此不断分解，蝶形结随着分解的次数成倍增加。每次迭代的蝶形类型比上一次蝶代增加一倍，数据点间隔也增大一倍。分解到最后，DFT只剩下运算只剩下了对蝶形结的运算，从而大大减小了运算量。

   频域抽取与时域抽取相同，同样是先奇偶分组，蝶形结构和时间抽取不一样但是蝶形个数一样，因此总的计算量也与时域抽取相同。

1.4Butterworth滤波器设计

滤波器可以使用Butterworth滤波器。Butterworth滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少,趋向负无穷大。

N阶Butterworth滤波器的衰减率为每倍频6 N分贝。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。滤波器阶数越高，在阻频带振幅衰减速度越快。此次程序设计中Butterworth滤波器的阶数为4。

二、系统框图：

图2.1 系统框图

三、代码实现和运行结果

3.1 导入数据

导入HWDDATE.M文件，为便于后续使用，将导入的变量s1重命名为HW。

3.2绘制原始信号时域波形

为了便于观察图像，将Plot窗口分成了4x1共四份

运行程序得到图3.1所示的波形。

图3.1 原始信号的时域波形

此时单从时域波形无法区分噪音和目标信号，因此需要对信号进行进一步处理。

3.3绘制原始信号的频域波形

先对原始信号进行FFT变换。再重新计算坐标数据。

随后绘制频域波形。     

运行程序得到图3.1所示的波形。

图3.2绘制原始信号的频域波形

经过观察，目标信号是94HZ到96HZ处的信号。

3.4设计Butterworth带通滤波器

先设计Butterworth带通滤波器，通过butter函数计算设计滤波器所需的参数，再用相应的参数通过filtfilt函数设计Butterworth滤波器。由图3.2知，目标信号是94HZ到96HZ处的信号因此设计滤波器阶数为4，滤波器的频率为94HZ到96HZ。

3.5绘制目标信号的时域波形

先对原始信号进行滤波，然后绘制目标信号的时域波形

绘制出的图像如图3.3

图3.3 目标信号的时域波形

经过观察，目标信号正常。

3.5绘制目标信号的频域波形

先对目标信号进行FFT变换，再绘制目标信号的频域波形。

这里仍使用到了fft。绘制出的图像如图3.4

图3.4 目标信号的频域波形

四、总结

代码总共绘制了四个图像，运行结果如下图：

图4.1 程序运行结果

对于该程序设计，最关键的是绘制频域波形，观察目标信号的频率。对于滤波器设计，MATLAB上有许多方法可以用来设计滤波器，比如filtfilt ( )函数butter ( )函数、designfilt( )函数等，每个函数都有自己的用法，有自己的参数，根据参数不同来设计不同的滤波器。但是直接写这些方法很不方便，我们可以借用MATLAB上的工具来设计滤波器，这里不再赘述。我设计的滤波器用到了filtfilt ( )函数,首先需要使用butter ( )函数计算滤波器参数，再使用butter ( )函数计算出的滤波器参数通过filtfilt ( ) 函数设计滤波器。就能得到目标信号了。