数字信号处理的MATLAB实现——快速傅里叶变换
文章目录
- 一、FFT函数
- 1.调用格式：
- 二、IFFT函数
- 1.调用格式
- 三、计算初始相位角与初始幅度
- 1.理论基础
- 2.调用格式
- 3.范例
- 四、频谱图中频率刻度的设置（以下讨论只考虑正频率）
- 1.N为奇数时
- 2.N为偶数时
- 3.范例
- 五、认识单频正弦信号相位谱
- 1.代码案例
- 2.分析
- 3.解决方法
- 六、消除趋势项
- 1.代码分析
- 2.问题分析
- 3.解决方法
- 七.将频谱图的纵坐标设置为分贝
- 1.方法：
- （1）y轴用对数坐标
- （2）转换为分贝值后再作图
- 八、频谱的栅栏现象
- 九、频谱的泄露现象
- 1.定义
- 2.通过加窗减小频谱泄漏
- （1）如何选取窗函数
- （2）N点DFT点数的选择
- （3）理想窗函数
- （4）解决加窗之后的幅值的衰减
- 十、分辨率
- 1.定义
- 十一、如何选择采样频率和信号长度
- 1.采样点数的选择
- （1）示例：
- （2）代码分析
- 十一.FFT中的补零问题
- 1.原理：
- 2.代码分析
- 3.注意：
- 十二、整周期采样
一、FFT函数
1.调用格式：
X=fft（x，N）
（1）当x长度大于N时，以x序列的长度做FFT
（2）当x长度小于N时，在x后补零至N长再做FFT
（3）当没有设置N时，以x序列的长度做FFT
二、IFFT函数
1.调用格式
x=ifft（X）
三、计算初始相位角与初始幅度
1.理论基础
当f0与FFT后频谱上的某根谱线相重合时，可求出对应的信号赋值和初始相角。但如果f0在两条谱线之间，就不能用这个方法来计算，这就是栅栏效应。
2.调用格式
（1）计算初始相位角
theta=abs（X（k））
其中X（k）为第k点的FFT值
（2）计算初始幅度
A=2*（abs（X（k）））/N
其中乘2是因为变换为单边谱，N为FFT点数，求出来的A为归一化幅度值。当频率在两条FFT谱线之间时，测量的幅度值有一定的误差。
3.范例
clear all; clc; close all;fs=1000; % 采样频率
N=1000; % 信号长度
t=(0:N-1)/fs; % 设置时间序列
f1=50; f2=65.75; % 两信号频率
x=cos(2*pi*f1*t)+cos(2*pi*f2*t); % 设置信号
X=fft(x); % FFT
Y=abs(X)*2/1000; % 计算幅值
freq=fs*(0:N/2)/1000; % 设置频率刻度
[A1, k1]=max(Y(45:65)); % 寻求第1个信号的幅值
k1=k1+44; % 修正索引号
[A2, k2]=max(Y(60:70)); % 寻求第1个信号的幅值
k2=k2+59; % 修正索引号
Theta1=angle(X(k1));
Theta2=angle(X(k2));
% 显示频率、幅值和初始相角
fprintf('f1=%5.2f A1=%5.4f Theta1=%5.4f\n',freq(k1),A1,Theta1);
fprintf('f2=%5.2f A2=%5.4f Theta2=%5.4f\n',freq(k2),A2,Theta2);% 作图
subplot 211; plot(freq,Y(1:N/2+1),'k'); xlim([0 150]);
xlabel('频率/Hz'); ylabel('幅值'); title('频谱图');
subplot 223; stem(freq,Y(1:N/2+1),'k'); xlim([40 60]);
xlabel('频率/Hz'); ylabel('幅值'); title('50Hz分量');
subplot 224; stem(freq,Y(1:N/2+1),'k'); xlim([55 75]);
xlabel('频率/Hz'); ylabel('幅值'); title('65.75Hz分量');四、频谱图中频率刻度的设置（以下讨论只考虑正频率）
1.N为奇数时
①freq=（0：(N-1)/2*fs/N）
③freq=linspace（0，fs/2-fs/N/2，(N-1)/2）
2.N为偶数时
①freq=（0：N/2-1）*fs/N
②freq=（0：N/2）*fs/N
③freq=linspace（0，fs/2，N/2+1）
3.范例
clear all; clc; close all;
fs=128; % 采样频率
N=128; % 信号长度
t=(0:N-1)/fs; % 时间序列
y=cos(2*pi*30*t); % 余弦信号
y=fft(y,N); % FFT
%freq=(0:N/2-1)*fs/N; % 按式(2-2-6c)设置正频率刻度
freq=linspace(0,fs/2,N/2)
% 作图
stem(freq,abs(y(1:N/2)),'k')
xlabel('频率(Hz)'); ylabel('幅值');
title('(b) 只有正频率刻度')
xlim([25 35]);
set(gcf,'color','w');
五、认识单频正弦信号相位谱
1.代码案例
%
% pr2_2_4
clear all; clc; close all;
fs=1000; % 采样频率
f0=50; % 信号频率
A=1; % 信号幅值
theta0=pi/3; % 信号初始相角
N=1000; % 信号长度
t=(0:N-1)/fs; % 设置时间序列
x=A*cos(2*pi*f0*t+theta0); % 设置信号
X=fft(x); % FFT
n2=1:N/2+1; % 设置索引号序列
freq=(n2-1)*fs/N; % 设置频率刻度
% 第一部分
THETA=angle(X(n2)); % 计算初始相角
Am=abs(X(n2)); % 计算幅值
ph0=THETA(51); % 计算信号的初始相角
fprintf('ph0=%5.4e %5.4e %5.4e\n',real(X(51)),imag(X(51)),ph0);
% 作图
%subplot 211; plot(freq,abs(X(n2))*2/N,'k');
%xlabel('频率/Hz'); ylabel('幅值')
%title('幅值谱图')
%subplot 212; plot(freq,THETA,'k')
%xlabel('频率/Hz'); ylabel('初始角/弧度')
%title('相位谱图')
%set(gcf,'color','w');
% 第二部分
Th=0.1; % 设置阈值
thetadex=find(Am<Th); % 寻找小于阈值的那线谱线的索引
THETA1=THETA; % 初始化THETA1
THETA1(thetadex)=0; % 对于小于阈值的那线谱线初始相位都为0
% 作图
figure(3)
%pos = get(gcf,'Position');
%set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-160)]);
plot(freq,THETA1,'k')
xlabel('频率/Hz'); ylabel('初始角/弧度')
title('相位谱图')
set(gcf,'color','w');
2.分析
因为每个频率都有初始相位，所以相位谱非常乱。
3.解决方法
可以通过设置有用频率分量最小幅值的方法，筛选出无用频率分量，将其赋值置为0。得到的频谱图如下图所示：
六、消除趋势项
1.代码分析
clear all; clc; close all;
load qldata.mat % 读入数据
N=length(y); % 数据长度
time=(0:N-1)/fs; % 时间刻度
% 第一部分
Y=fft(y); % FFT
n2=1:N/2+1; % 取正频率索引序列
freq=(n2-1)*fs/N; % 频率刻度
% 作图
subplot 211; plot(time,y,'k'); ylim([0 15]); grid;
title('有趋势项的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(Y(n2)),'k')
title('有趋势项的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');
% 第二部分
x=detrend(y); % 消除趋势项
X=fft(x); % FFT
% 作图
figure
subplot 211; plot(time,x,'k'); ylim([-5 5]); grid;
title('消除趋势项后的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(X(n2)),'k');
title('消除趋势项后的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');
2.问题分析
可以看到FFT频谱图中什么都没有，几乎都为0。这是因为波形中有一个明显的趋势项。正是因为该趋势项的存在使得频谱分析中有很大的直流分量。
3.解决方法
为了消除趋势项，可采用detrend函数：
xd=detrend(x)
消除趋势项后的FFT，如图所示：
七.将频谱图的纵坐标设置为分贝
1.方法：
（1）y轴用对数坐标
调用semilog函数
semilogx(x,y) %x轴为对数坐标
semilogy(x,y) %y轴为对数坐标
loglog(x,y) %x，y轴都为对数坐标
（2）转换为分贝值后再作图
令Y=20*log10(abs(Y))
八、频谱的栅栏现象
九、频谱的泄露现象
1.定义
对于连续信号的采样序列进行DFT运算，由于时间长度取有限值，即将信号截断，使信号的带宽被扩展了，这种现象称为泄露。
2.通过加窗减小频谱泄漏
（1）如何选取窗函数
选择第一旁瓣衰减大，旁瓣峰值衰减快的窗函数有利于缓解截断过程中产生的频谱泄漏问题。
（2）N点DFT点数的选择
其中K为窗函数的主瓣宽度与矩形窗的主瓣宽度之比。
（3）理想窗函数
①窗函数的频谱尽可能窄，以提高谱估计时的频域分辨率和减小泄漏。
②尽量减小窗函数频谱的最大旁瓣的相对幅度，以使旁瓣高度随频率尽快衰减。
（4）解决加窗之后的幅值的衰减
为了修正加窗造成的幅值衰减，可在加窗FFT后的频谱进行修成，乘上窗函数的恢复系数
十、分辨率
1.定义
（1）物理分辨率：谱分析中将信号中两个靠的很近的谱峰仍能保持分辨的能力。
（2）计算分辨率使用DFT时，在频率轴上的所能得到的最小频率间隔。
十一、如何选择采样频率和信号长度
1.采样点数的选择
（1）示例：
（2）代码分析
①代码块：
clear all; clc; close all;
M=256; fs=10; % 设置数据长度M和采样频率fs
f1=1; f2=2.5; f3=3; % 设置3个正弦信号的频率
t=(0:M-1)/fs; % 设置时间序列
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); % 计算出信号波形X1=fft(x,20); % FFT变换
X2=fft(x,40);
X3=fft(x,128);
freq1=(0:10)*fs/20; % 计算3个信号在频域的频率刻度
freq2=(0:20)*fs/40;
freq3=(0:64)*fs/128;
% 作图
plot(freq1,abs(X1(1:11)),'k',freq2,abs(X2(1:21)),'k',freq3,abs(X3(1:65)),'k');
%legend('N=20','N=40','N=128');
title('不同N值的DFT变换');
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');
②结果：N的点数不同，正弦信号的幅值也不同
十一.FFT中的补零问题
1.原理：
在FFT时经常会对输入序列补零，最常见的是补零后把N设置为2的k次幂。补零实际上是对FFT后的X(k)进行插值，补零后能改善栅栏效应，使原先不清晰的谱线显现出来。
2.代码分析
%
% pr2_2_12
clear all; clc; close all;fs=200; % 采样频率
f1=30; f2=65.5; % 两信号频率
N=200; % 信号长度
n=1:N; % 样点索引
t=(n-1)/fs; % 时间刻度
x=cos(2*pi*f1*t)+cos(2*pi*f2*t); % 信号X1=fft(x); % 按N点进行FFT
freq1=(0:N/2)*fs/N; % N点时正频率刻度
X1_abs=abs(X1(1:N/2+1))*2/N; % 信号幅值L=2*N; % 补零后FFT长度
X2=fft(x,L); % 按L长进行FFT
freq2=(0:L/2)*fs/L; % L点时频率刻度
X2_abs=abs(X2(1:L/2+1))*2/N; % 信号幅值
% 作图
subplot 211; plot(freq1,X1_abs,'k');
grid; ylim([0 1.2]);
xlabel('频率/Hz'); ylabel('幅值');
title('(a) 补零前FFT谱图')
subplot 212; plot(freq2,X2_abs,'k');
grid; ylim([0 1.2]);
xlabel('频率/Hz'); ylabel('幅值');
title('(b) 补零后FFT谱图')
set(gcf,'color','w');
由于在FFT上内插，所以可以看到主峰的周围有波纹震荡，这是由于泄露造成的。实际上补零前的信号并不是没有泄漏，而是在谱线的各个频点上正好与谱线的零点相重合。
3.注意：
补零不能增加物理分辨率，只能增加计算分辨率。要增加物理分辨率，唯一的方法是增加信号的有效长度。
十二、整周期采样
数字信号处理的MATLAB实现——快速傅里叶变换相关推荐
- dft对称性 matlab实验,数字信号处理实验 matlab版 离散傅里叶变换的性质
数字处理实验 matlab版 山大学生最适用 本人自己写的 因为时间比较久了 不能完全保证出现代码都能运行 但95%还是能保证的 谢谢 实验13 离散傅里叶变换的性质 (完美格式版,本人自己完成,所有 ...
- 数字信号处理实验matlab版答案刘舒帆,数字信号处理实验（MATLAB版） 刘舒帆,费诺,陆辉 西安电子科技大学出版社 9787560620060...
商品描述: 基本信息 书名:数字信号处理实验(MATLAB版) 原价:31.00元 作者:刘舒帆,费诺,陆辉 著 出版社:西安电子科技大学出版社 出版日期:2013-7-1 ISBN:97875606 ...
- matlab的dft谱分析,数字信号处理基于matlab(用DFT作谱分析,窗函数的设计)
数字信号处理基于matlab(用DFT作谱分析,窗函数的设计) 1实验一用DFT作谱分析X11111X212344321N108X3COSN1PI/4N208X4SINN2PI/8FIGURESUBP ...
- 【数字信号处理及MATLAB实践】
数字信号处理及MATLAB实践 第一章 信号.连续时间周期信号的傅里叶级数和频谱分析 文章目录 数字信号处理及MATLAB实践 前言 1.1 信号的时域分析-波形的产生和信号的基本运算及MATLAB实 ...
- 全相位数字信号处理方法及matlab实现,数字信号处理及matlab实现_实验报告册.doc...
数字信号处理及matlab实现_实验报告册.doc 数字信号处理及MATLAB实现实验报告实验人孙敬贤实验1离散时间信号产生及频谱分析一.实验目的㈠掌握MATLAB产生常用离散时间信号的产生方法.㈡掌 ...
- matlab fftshift_数字信号处理没有Matlab？用Python一样很爽
通常,在数字信号处理时,我们避不开matlab这个工具,因其它的强大的功能受到广大工程师的好评,也一直都是业界的不二之选.但是,matlab毕竟是商业软件,公司里如果使用的话,就需要支付高昂的费用.即 ...
- matlab抽样仿真混叠图,数字信号处理及MATLAB仿真__前言
前言 历史背景 像许多电子工程课程一样,数字信号处理(Digital Signal Processing,DSP)最初是一门研究生课程,近30年来,其逐渐向本科课程渗透,成为电子与计算机工程的本科课程 ...
- matlab如何进行数字信号处理,使用MATLAB进行数字信号处理-第2部分
这篇文章来源于DevicePlus.com英语网站的翻译稿. 在Arduino DSP系列的第二部分中,我们将继续深入研究数字信号处理的基础知识.我们将学习数字滤波器的特性以及在MATLAB中处理信号 ...
- Python学习-Scipy库信号处理signal(过滤、快速傅里叶变换、信号窗函数、卷积)
Python学习-Scipy库信号处理signal 目录 1.过滤:以某种方式修改输入信号 2.快速傅里叶变换 3.信号窗函数 4.卷积 导入库 import matplotlib.pyplot as ...
最新文章
- 量子位MEET大会报名开启！各领域头部玩家集结，AI年度榜单揭晓，在这里预见智能科技新未来...
- 已知服务器ftp的账号密码，求解数据库表的内容
- 联想天工 802.1x认证 主程序
- bzoj 2563 贪心 思想
- opencv中匹配点对的坐标提取
- PHPCMS之 列表和内容页
- matlab 四种取整函数(fix floor ceil round)的区别
- php laravel 理解,程序员-说一下PHP框架Laravel，如何理解她的思想
- Python（八）：条件与循环
- 32 道常见的 Kafka 面试题
- windows系统——更改系统关机音效
- 【隐形的翅膀】基于钉钉工作流的人事评价信息采集案例(2)：钉钉智能表单、OA审批、自动任务功能对比
- css 大于号 标签_css选择器 ~ （波浪号）、+（加号）、（大于号）的用法解析和举例...
- At least one JAR was scanned for TLDs yet contained no TLDs.问题解决方式
- 愤怒的小鸟 高清完整版下载
- 【程序源代码】小程序商城系统（CoreShop）
- fullCalendar改造计划之带农历节气节假日的万年历
- Feign - 独立使用 - 替代HttpClient
- 微软enchange服务器安装,安装 Exchange Server 2010
- css卷轴动画小程序,微信小程序动画两种实现方式