MATLAB中功率谱密度计算pwelch函数使用详解

MATLAB中功率谱密度计算pwelch函数使用详解

目录

前言

一、pwelch函数简介

二、pwelch函数参数说明

三、pxx = pwelch(x)示例

四、[pxx,f]=pwelch(x,window,noverlap,nfft,fs)示例

四、[pxx,f] = pwelch(x,window,noverlap,nfft,fs,freqrange,spectrumtype)示例

五、多通道功率谱估计

六、参考资料

总结


前言

        详细介绍MATLAB中功率谱密度计算pwelch函数的使用方法,介绍如何使用该函数及输入各个参数的含义,手把手用代码教你学习pwelch函数,文中附有代码,足够pwelch函数入门了。


提示:以下是本篇文章正文内容,希望能帮助到各位,转载请附上链接。

一、pwelch函数简介

         MATLAB中的pwelch函数是一种用于快速估计信号功率谱密度的工具,也可以计算信号的功率谱,通过阅读该函数使用说明会发现功率谱和功率谱密度是两个不同的概念,要注意一下,在很多教材上都称功率谱和功率谱密度是同一个概念,这是错的,不要被误导。

        pwelch函数可以只对一个信号进行功率谱密度估计,也可以同时对多个信号进行功率谱密度估计,简而言之,就是其输入信号那个参数可以是向量,也可以是矩阵。

二、pwelch函数参数说明

        函数说明文档里面pwelch函数调用的句话格式很多,我们不用关心那么多,关心如下几个格式和默认参数是怎么一回事就可以了。

pxx = pwelch(x)

[ pxx, f ] = pwelch(x,window,noverlap,nfft,fs)

[ pxx, f ] = pwelch(x,window,noverlap,nfft,fs,freqrange)

[ pxx, f ] = pwelch(x,window,noverlap,nfft,fs,freqrange,spectrumtype)

x--------输入信号,指定为行向量或列向量,或矩阵。 如果 x 是矩阵,则其被视为独立通道

Window--------窗口,指定为行向量或列向量或整数。 如果 window 是一个向量,pwelch 将 x 划分为长度等于 window 长度的重叠段,然后将每个信号段与 window 中指定的向量相乘。 如果window是整数,则将pwelch分成长度等于整数值的段,并使用等长的汉明窗。 如果x的长度不能精确地划分为具有noverlap数量的重叠样本的整数段,则x被相应地截断。 如果将 window 指定为空,则使用默认的 Hamming 窗来获取 x 的 8 段,其中具有 noverlap 重叠样本。

noverlap----------重叠样本的数量,指定为小于窗口长度的正整数。 如果省略 noverlap 或 noverlap 指定为空,则使用一个值来获得段之间 50% 的重叠,即默认50%重叠

nfft--------DFT 点数,指定为正整数。 对于实值输入信号 x(PSD 估计值),如果 nfft 为偶数,则 pxx 的长度为 (nfft/2 + 1);如果 nfft 为奇数,则 pxx 的长度为 (nfft + 1)/2。 对于复值输入信号 x,PSD 估计的长度始终为 nfft。 如果 nfft 指定为空,则使用默认的 nfft。如果 nfft 大于Window长度,则数据用零填充。 如果 nfft 小于Window长度,则使用 datawrap包装该段,使长度等于nfft。建议两者相等。

fs-------采样频率,指定为正标量。 采样率是单位时间内的采样数。 如果时间单位是秒,那么采样率的单位是Hz。

freqrange--------------PSD 估计的频率范围,指定为“单边”、“双边”或“中心”之一。 对于实值信号,默认值为“单侧”;对于复值信号,默认值为“双侧”。 每个选项对应的频率范围为

        'oneside' — 返回实值输入信号 x 的单侧 PSD 估计。 如果 nfft 为偶数,则 pxx 的长度为 nfft/2 + 1, 如果 nfft 为奇数,则 pxx 的长度为 (nfft + 1)/2,间隔为 [0,π) rad/sample。 当 fs 可选指定时,对于偶数和奇数长度 nfft,相应的间隔分别为 [0,fs/2] 周期/单位时间和 [0,fs/2) 周期/单位时间。该函数将除 0 和奈奎斯特频率之外的所有频率的功率乘以 2,以保持总功率不变。

        'twoside' - 返回实值或复值输入 x 的两侧 PSD 估计值。 在这种情况下,pxx 的长度为 nfft,并在区间 [0,2π) rad/sample 上计算。 当指定fs时,间隔为[0,fs)周期/单位时间。

        'centered' - 返回实值或复值输入 x 的中心两侧 PSD 估计值。 在这种情况下,pxx 的长度为 nfft,并且在偶数长度 nfft 的间隔 (–π,π] rad/sample 和奇数长度 nfft 的 (–π,π) rad/sample 区间内计算。当 fs 可选地指定时,相应的对于偶数和奇数长度 nfft,间隔分别为 (–fs/2, fs/2] 周期/单位时间和 (–fs/2, fs/2) 周期/单位时间。

spectrumtype------------功率谱类型,指定为“psd”或“power”之一。 默认“psd”,将返回功率谱密度。 指定“power”可通过窗口的等效噪声带宽来缩放 PSD 的每个估计值。 使用“power”选项可以获得每个频率的功率估计值。

三、pxx = pwelch(x)示例

        创建一个角频率为π/4 rad的正弦波,外加N(0,1)白噪声。信号的长度是320个样本。得到其Welch PSD估计。

clc;
clear;
close all;rng defaultn = 0:319;
x = cos(pi/4*n)+randn(size(n));pxx = pwelch(x);

可见,默认横轴是归一化的角频率,纵轴是取了10log10( )的dB功率谱密度。

关于归一化频率,参考:滤波器设计中的频率归一化问题_滤波器归一化频率-CSDN博客

解释如下:

        信号处理工具箱中经常使用的频率是Nyquist频率,它被定义为采样频率的一半,在滤波器的结束选择和设计当中的截止频率均使用Nyquist频率进行归一化处理。

   例如,对于一个采样频率为1000Hz的系统,300Hz的归一化即为300/500=0.6。归一化频率的范围在[0,1]之间。如果要将归一化频率转换为角频率,则将归一化频率乘以pi;如果将归一化频率转换成Hz,则将归一化频率乘以采样频率的一半。

        采样率的一半是最高频率,认为是1,那么真实频率和最高频率的比值就是归一化频率!也叫数字频率。
     将信号分成长度为nsc=⌊Nx/4.5⌋。这个动作相当于将信号分成尽可能长的段,以获得接近但不超过8个重叠50%的段。使用汉明窗口显示各部分。指定相邻部分之间50%的重叠

要计算FFT,使用max(256,2^p),其中p=[log2nsc⌉。

Nx = length(x);
nsc = floor(Nx/4.5);
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));t = pwelch(x,hamming(nsc),nov,nff);maxerr = max(abs(abs(t(:))-abs(pxx(:))))

默认分成8段,每段的长度为Nx/4.5=320/4.5=71.1111,舍去多余的数据,分段长度为71,并用同等长度的汉明窗;重叠长度为50%,则nov=floor(71/2)=35;DFT的点数取每段长度最接近的2的整数次幂和256的最大值,最接近的2的整数次幂是比每段长度长的最接近的2的整数次幂,所以DFT计算的时候如果DFT点数大于每段长度。会自动补0。

对于实值信号,默认值为“单侧”PSD,所以计算DFT点数为256,估计的PSD长度只有256/2+1=129点的长度。

四、[pxx,f]=pwelch(x,window,noverlap,nfft,fs)示例

        创建一个由100Hz正弦波和加性N(0,1)白噪声组成的信号。采样率为1khz,信号持续时间为5秒。使用500个样本和300个重叠样本的段长度,使用500个DFT点。

clc;
clear;
close all;rng defaultfs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));[pxx,f] = pwelch(x,500,300,500,fs);plot(f,10*log10(pxx))xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

sum(abs(x).^2)/length(x) %信号功率
P = fs/2/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分

验证功率,对功率谱密度积分发现和信号功率几乎相等。

增大点数

clc;
clear;
close all;rng defaultfs = 10000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));[pxx,f] = pwelch(x,5000,3000,5000,fs);plot(f,10*log10(pxx))xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')sum(abs(x).^2)/length(x) %信号功率
P = fs/2/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分

增大段长度,会发现功率谱估计的准,功率谱密度积分的值更接近信号的功率。

四、[pxx,f] = pwelch(x,window,noverlap,nfft,fs,freqrange,spectrumtype)示例

         创建一个由100Hz正弦波和加性N(0,1)白噪声组成的信号。采样率为1khz,信号持续时间为5秒。使用500个样本和300个重叠样本的段长度,使用500个DFT点。

clc;
clear;
close all;rng defaultfs = 1000;
t = 0:1/fs:10-1/fs;x = cos(2*pi*100*t)+randn(size(t))+1;[pxx,f] = pwelch(x,500,300,500,fs,'centered');plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
E=sum(abs(x).^2)  %信号能量
P=sum(abs(x).^2)/length(x) %信号功率
P = fs/length(pxx) * (sum(pxx) - 0.5*(pxx(1) + pxx(end))) %梯形法积分�[pxx,f] = pwelch(x,500,300,500,fs,'centered','power');
figure(2)
plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

使用参数centered得到双边功率谱密度:

使用参数power得到双边功率谱(不是双边功率谱密度):

观察图可知,与信号设置的功率吻合。

取对数画出如下:

        可见,将功率谱密度积分和信号功率相等;观察功率谱图可知,每个频率对应的点显示了该频点的功率大小。

五、多通道功率谱估计

        在加性N(0,1)高斯白噪声中生成由三个正弦波组成的多通道信号的1024个样本。正弦波的频率分别为100、200、300Hz,采样频率为1000Hz。使用Welch的方法估计信号的PSD并绘制它。

clc;
clear;
close all;rng defaultfs = 1000;
t = 0:1/fs:5-1/fs;
f = [100;200;300];
x = cos(2*pi*f*t)' + randn(length(t),3);[pxx,f] = pwelch(x,500,300,500,fs);plot(f,10*log10(pxx));xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

六、参考资料

Welch’s power spectral density estimate - MATLAB pwelch- MathWorks 中国

https://download.csdn.net/download/m0_66360845/89084990icon-default.png?t=N7T8https://download.csdn.net/download/m0_66360845/89084990


总结

        以上就是今天要讲的内容,本文仅仅简单介绍了功率谱密度(功率谱)绘制函数pwelch函数的使用,希望对各位小伙伴有所帮助。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/bicheng/6294.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

Springboot(SSM)项目实现数据脱敏

目录 一、引入hutool的依赖 二、sql脚本 三、自定义注解代码 3.1 自定义注解 3.2 自定义一个枚举,用于定义脱敏的类型 3.3 序列化 四、使用脱敏注解 4.1 Person.java 4.2 controller 4.3 dao 五、源代码参考 一、引入hutool的依赖 <dependency><groupId>…

Linux下Zip命令无法压缩软链接指向的实际文件的解决方案

Linux下Zip命令无法压缩软链接指向的实际文件的解决方案 一、前言 最近在作横向课题&#xff0c;在进行RKNN1808的目标检测C推理环境部署时&#xff0c;遇到了视频和图像的处理问题。出现这些问题&#xff0c;最好是找谁呢&#xff1f;找OpenCV。 但是本身开发板是个空架子&a…

FIFO Generate IP核使用——Native读写接口信号详解

Native FIFO接口信号是用于FIFO IP核与外部电路进行通信的信号。当FIFO支持独立的写和读时钟时&#xff0c;这些信号可以包括标准端口和可选端口。 1 当FIFO具有独立时钟时的接口信号 当FIFO具有独立的时钟时&#xff0c;其接口信号会相应地有所变化。特别是关于复位信号rst…

计算机等级考试2级(Python)知识点整理

计算机等级考试2级&#xff08;Python&#xff09;知识点整理 1.基础知识点&#xff08;记忆、理解&#xff09; 第1讲Python概述 01. 源代码 02. 目标代码 03. 编译和解释 04. 程序的基本编写方法 第2讲 Python语言基础&#xff08;一&#xff09; 01. 用缩进表示代码…

[数据结构]——非比较排序—计数排序

该篇文章 所涉及代码收录仓库&#xff1a;登录 - Gitee.com 目录 1.非比较排序——计数排序 2.最终实现 1.解析 2.以int a[] { 1,3,9,1,5,1,2,3,-5,-5,-2 };为例&#xff0c;手撕分析 3.代码实现 4.计数排序具有以下主要特性&#xff1a; 1.非比较排序——计数排序 思想…

Linux IP Forwarding路由转发实验

linux 路由转发功能 Linux 操作系统具备路由转发功能&#xff0c;路由功能是指 Linux 操作系统提供的路由管理和转发功能&#xff0c;它允许 Linux 主机在网络中正确地转发数据包&#xff0c;并确保数据包能够达到其目的地。 出于安全考虑&#xff0c;Linux系统默认是禁止数据…

ICode国际青少年编程竞赛- Python-1级训练场-for循环入门

ICode国际青少年编程竞赛- Python-1级训练场-for循环入门 1、 for i in range(4):Dev.step(4)Dev.turnLeft()2、 for i in range(3):Dev.step(6)Dev.turnRight()3、 for i in range(3):Dev.turnRight()Dev.step(2)Dev.turnLeft()Dev.step(-3)4、 for i in range(4):Dev…

附录E:董事会

< 回到目录 附录E&#xff1a;董事会 当董事会和CEO配合默契时&#xff0c;董事会成员扮演着CEO副手的角色。董事会还对CEO进行监督&#xff0c;如果CEO出现失职的情况&#xff0c;董事会将解雇他&#xff0c;并另请高明。这是董事会的职责所在。很遗憾&#xff0c;大多数董…

谷神前端组件增强:二级列表

createGuthonOptBtnsElement /** * 根据提供的row、options为子列表创建并返回谷神操作按钮元素。 * * param {object} row - agGrid单元格渲染函数参数。* param {array} options - 谷神操作按钮配置项数组。* returns {Element} - 返回创建的DOM元素。 */ function cr…

数字旅游以科技创新为核心竞争力:推动旅游服务的智能化、高效化,满足游客日益增长的旅游需求

一、引言 随着科技的飞速发展&#xff0c;数字旅游作为旅游业与信息技术结合的产物&#xff0c;正以其独特的魅力改变着传统旅游业的格局。科技创新作为数字旅游的核心竞争力&#xff0c;不仅推动了旅游服务的智能化、高效化&#xff0c;更满足了游客日益增长的旅游需求。本文…

香港理工大学内地事务总监陆海天教授确认出席“边缘智能2024 - AI开发者峰会”并发表主题演讲

隨著AI技術的日新月異&#xff0c;我們正步入一個邊緣計算智能化與分布式AI相互融合的新紀元。這一變革不僅推動了分布式智能創新應用的飛速發展&#xff0c;還使得邊緣智能——這一結合邊緣計算和智能技術的新興領域&#xff0c;逐漸成為引領AI發展的重要力量。通過其分布式和…

clang:在 Win10 上编译 MIDI 音乐程序(二)

先从 Microsoft C Build Tools - Visual Studio 下载 1.73GB 安装 "Microsoft C Build Tools“ 访问 Swift.org - Download Swift 找到 Windows 10&#xff1a;x86_64 下载 swift-5.10-RELEASE-windows10.exe 大约490MB 建议安装在 D:\Swift\ &#xff0c;安装后大约占…

llmperf测试大模型API性能

llmperf测试大模型API性能 llmperf是一个用来评估LLM API性能的工具。 官方仓库地址&#xff1a;https://github.com/ray-project/llmperf 1. 安装准备 脚本依赖python3环境&#xff0c;测试前客户端安装python3&#xff0c;本文使用python版本为3.8。 # 创建一个python虚拟环境…

走进香港美食宛如走进香港电影

&#xff08;1&#xff09; 过去蔡澜有个节目&#xff0c;专门介绍香港美食&#xff0c;身边美女相伴、眼里美景相随。 过去离香港海关近&#xff0c;有时候散步都能走到那里&#xff0c;打车时车都不蹦字儿。那时候精神头儿真好&#xff0c;周六一早6点就起来拖着大箱子过关&a…

React 第十八章 Hook useImperativeHandle

React 的 useImperativeHandle 是一个自定义 Hook。该 Hook 一般配合 React.forwardRef 使用&#xff0c;主要用于在父组件中可以使用 ref 获取子组件暴露出来的属性和方法。 useImperativeHandle 接受两个参数&#xff1a;ref 和创建子组件属性和方法的函数。 function Chil…

人机协同中的分布式中心化态势感知

一、集中式和分布式 集中式和分布式是两种不同的系统结构和管理方式。集中式系统是指所有计算机资源和数据都集中在一个中心服务器或主机上&#xff0c;所有的计算和数据处理都由该中心服务器来完成。而分布式系统是指计算机资源和数据分布在不同的计算机节点上&#xff0c;通…

软件测试(实验五)——Jmeter的使用

目录 实验目的 一、使用JMeter演示取样器、监听器、配置元件、断言的使用&#xff1b; 1、取样器 2、监听器 3、配置元件的使用 ① 用户定义的变量 ②HTTP信息头管理器 ③HTTP请求默认值 ④CSV数据文件设置 4、断言 ①响应断言 ②JSON断言 ③断言持续时间 二、使用…

「 网络安全常用术语解读 」SBOM主流格式SWID详解

国际标准化组织&#xff08;ISO&#xff09;和国际电工委员会&#xff08;International Electrotechnical Commission&#xff0c;IEC&#xff09;发布了ISO/IEC 19770-2软件标识&#xff08;Software Identification&#xff0c;SWID&#xff09;标签标准&#xff0c;该标准定…

【LAMMPS学习】八、基础知识(5.6)绝热核/壳模型

8. 基础知识 此部分描述了如何使用 LAMMPS 为用户和开发人员执行各种任务。术语表页面还列出了 MD 术语&#xff0c;以及相应 LAMMPS 手册页的链接。 LAMMPS 源代码分发的 examples 目录中包含的示例输入脚本以及示例脚本页面上突出显示的示例输入脚本还展示了如何设置和运行各…

4.2 JavaScript语法

4.2.1 JavaScript大小写 在JavaScript中大小写是严格区分的&#xff0c;无论是变量、函数名称、运算符和其他语法都必须严格按照要求的大小写进行声明和使用。例如变量hello与变量HELLO会被认为是完全不同的内容。 4.2.2 JavaScript分号 很多编程语言&#xff08;例如C、Java和…