信号频谱分析、功率谱分析、倒谱分析、小波分析
PDF获取路径
1.频谱分析
在本科信号系统课程中学习过傅里叶变换,可将信号时域波形转化为频域。为什么要进行域转换呢?因为大部分信号在传输过程中可能会受到外界因素的干扰(可以理解为"噪声"),这种干扰在时域上表现得不太明显,因此可以通过傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号(信号的频谱)。
傅立叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号,以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。和傅立叶变换算法对应的是反傅立叶变换算法。该反变换从本质上说也是一种累加处理,这样就可以将单独改变的正弦波信号转换成一个信号。
clear; clc;
t_s = 0.01; % 采样周期, 采样频率fs为100Hz
t_start = 0.5; % 起始时间
t_end = 5; % 结束时间
t = t_start : t_s : t_end;
y = 1.5*sin(2*pi*5*t) + 3*sin(2*pi*20*t) + randn(1,length(t)); % 生成信号%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频谱%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_f = fft(y); % 傅里叶变换
figure;
subplot(5,1,1);
plot(t,y); title('original signal'); % 绘制原始信号图
Druation = t_end - t_start; % 计算采样时间
Sampling_points = Druation/t_s + 1; % 采样点数,fft后的点数就是这个数
f_s = 1/t_s; % 采样频率
f_x = 0:f_s/(Sampling_points -1):f_s; % 注意这里和横坐标频率对应上了,频率分辨率就是f_s/(Sampling_points -1)
t2 = f_x-f_s/2;
shift_f = abs(fftshift(y_f));
subplot(5,1,2);
plot(f_x,abs(y_f));title('fft transform');
subplot(5,1,3);
plot(f_x-f_s/2,shift_f);title('shift fft transform'); % 将0频率分量移到坐标中心
subplot(5,1,4);
len_t2 = floor(length(t2)/2);
len_f = floor(length(shift_f)/2);
plot(t2(len_t2:len_t2*2),shift_f(len_f:len_f*2));
title('shift fft transform'); % 保留正频率部分
subplot(5,1,5);
len_fx = floor(length(f_x)/2);
plot(f_x(1:len_fx),abs(y_f(1:len_fx)));
title('fft cut'); % 直接截取fft结果的前半部分
第一幅图为生成信号的时域波形图,频率成分为5、20Hz,含一定程度噪声干扰。
第二幅图为该信号经傅里叶变换得到的频谱图,中心轴对称;
第三幅图为图二频谱搬移情况,将0频率分量搬到中心;
第四幅图为图三正频率部分;
第五幅图为图二结果的前半部分。
2. 功率谱分析
另外,可以从含噪信号功率谱中观察各个子频率强度;求取功率谱有fft法、周期图法和自相关法三种方式。
% 1.(傅立叶变换的平方)/(区间长度);->直接法
% 2.自相关函数的傅里叶变换;->相关函数法
clear; clc;
Fs = 1000; % 采样频率
nfft = 1000; % 采样点数% 产生序列
n = 0:1/Fs:1;
xn = cos(2*pi*100*n) + 3*cos(2*pi*200*n) + (randn(size(n)));
figure;
subplot(5,1,1); plot(xn);
title('加噪信号'); xlim([0,1000]); grid on;
% FFT
Y = abs(fft(xn,nfft));
subplot(5,1,2); plot((10*log10(Y(1:nfft/2))));
title('FFT'); xlim([0,500]); grid on;
% 法1:FFT直接平方
Y2 = Y.^2/(nfft);
subplot(5,1,3); plot(10*log10(Y2(1:nfft/2)));
title('直接法'); xlim([0,500]); grid on;
% 法2:周期图法
window = boxcar(length(xn)); % 矩形窗
[psd1,f] = periodogram(xn,window,nfft,Fs);
psd1 = psd1 / max(psd1); % 归一化
subplot(5,1,4); plot(f,10*log10(psd1));
title('周期图法'); ylim([-60,10]); grid on;
% 法3:自相关结果
cxn = xcorr(xn,'unbiased'); % 计算自相关函数
% 法4:自相关法
CXk = fft(cxn,nfft);
psd2 = abs(CXk);
index = 0:round(nfft/2-1);
k = index*Fs/nfft;
psd2 = psd2/max(psd2); % 归一化
psd2 = 10*log10(psd2(index+1)); % 拉高
subplot(5,1,5); plot(k,psd2); title('间接法'); grid on;
% 取对数的目的是使那些振幅较低的成分相对高振幅成分得以拉高,
% 以便观察掩盖在低幅噪声中的信号特征。
从图中看出通过自相关法得到的功率谱成分(100Hz,200Hz)更加明显,其它频带干扰较少。
3. 倒谱分析
科学上网
blog
clear; clc;
sf = 1000; % 采样频率
nfft = 1000; % 采样点数
x = 0:1/sf:5;
y1 = 10*cos(2*pi*5*x)+7*cos(2*pi*10*x)+5*cos(2*pi*20*x)+0.5*randn(size(x));
y2 = 20*cos(2*pi*50*x)+15*cos(2*pi*100*x)+25*cos(2*pi*200*x)+0.5*randn(size(x));
for i = 1:length(x)y(i) = y1(i)*y2(i); % 被调制
end
figure;
subplot(3,3,1)
plot(y1);xlim([0 5000]);title('y1');
subplot(3,3,2)
plot(y2);xlim([0 5000]);title('y2');
subplot(3,3,3)
plot(y);xlim([0 5000]);title('y=y1*y2');t = 0:1/sf:(nfft-1)/sf;
nn = 1:nfft;subplot(3,3,4)
Y = abs(fft(y1,nfft)); % 1024点fft
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y_1');
ylabel('幅值');xlim([0 300]);
grid on;subplot(3,3,5)
Y = abs(fft(y2,nfft));
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y_2');
ylabel('幅值');xlim([0 300]);
grid on;subplot(3,3,6)
Y = abs(fft(y,nfft));
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y');
ylabel('幅值');xlim([0 300]);
grid on;subplot(3,3,7)
z = rceps(y);
plot(t(nn),abs(z(nn)));
title('z=rceps(y)');ylim([0 0.3]);
xlabel('时间(s)');
ylabel('幅值');
grid on;
subplot(3,3,8)
yy = real(ifft(log(abs(fft(y))))); % 信号→傅里叶→对数→傅里叶逆变换
plot(t(nn),abs(yy(nn)));
title('real(ifft(log(abs(fft(y)))))');ylim([0 0.3]);
xlabel('时间(s)');
ylabel('幅值');
grid on;
可以发现第7幅图跟第8幅图一样,说明了倒谱等价于信号→傅里叶→对数→傅里叶逆变换。
进阶
举一个例子:
通过倒谱达到语音分离效果;
注意:这里的时间并不是严格意义上的时间,而是伪频率。
load mtlb% soundsc(mtlb,Fs);
% To hear, type soundsc(mtlb,Fs)timelag = 0.23; % 时延
delta = round(Fs*timelag);
alpha = 0.5; % 衰减系数orig = [mtlb;zeros(delta,1)];
echo = [zeros(delta,1);mtlb]*alpha;mtEcho = orig + echo;
pause(0.5);soundsc(mtEcho,Fs);t = (0:length(mtEcho)-1)/Fs;
figure;
subplot(2,1,1)
plot(t,[orig echo])
legend('Original','Echo')subplot(2,1,2)
plot(t,mtEcho)
legend('Total')
xlabel('Time (s)')
c = rceps(mtEcho);[px,locs] = findpeaks(c,'Threshold',0.2,'MinPeakDistance',0.2);clf
plot(t,c,t(locs),px,'o')
xlabel('Time (s)')dl = locs(2)-1;mtNew = filter(1,[1 zeros(1,dl-1) alpha],mtEcho);subplot(2,1,1)
plot(t,orig)
legend('Original')subplot(2,1,2)
plot(t,mtNew)
legend('Filtered')
xlabel('Time (s)')
% To hear, type soundsc(mtNew,Fs)
soundsc(mtNew,Fs)
于是,两个语音段得到了分离;
4. 小波分析
clear; clc;
t_s = 0.005; % 采样周期
t_start = 0.001; % 起始时间
t_end = 10; % 结束时间
t = t_start : t_s : t_end;
y = 10*sin(2*pi*2*t)+3*sin(2*pi*10*t)+1*sin(2*pi*20*t)+3*randn(1,length(t)); %生成信号
len = length(y);
% 生成突变信号
y2 = 50*sin(2*pi*50*t); % 高振幅高频率 突变信号
for i = 1: lenif i>=601&&i<=604y(i) = y(i)+y2(i);elsey(i) = y(i);end
end
figure;
plot(y) % 绘制原始信号
[c,l] = wavedec(y,5,'db5');
% 重构1~5层细节函数
d5 = wrcoef('d',c,l,'db5',5);
d4 = wrcoef('d',c,l,'db5',4);
d3 = wrcoef('d',c,l,'db5',3);
d2 = wrcoef('d',c,l,'db5',2);
d1 = wrcoef('d',c,l,'db5',1);
% 重构1~5层近似函数
a5 = wrcoef('a',c,l,'db5',5);
a4 = wrcoef('a',c,l,'db5',4);
a3 = wrcoef('a',c,l,'db5',3);
a2 = wrcoef('a',c,l,'db5',2);
a1 = wrcoef('a',c,l,'db5',1);
figure
subplot(5,2,1);
plot(a1)
subplot(5,2,2);
plot(d1)
subplot(5,2,3);
plot(a2)
subplot(5,2,4);
plot(d2)
subplot(5,2,5);
plot(a3)
subplot(5,2,6);
plot(d3)
subplot(5,2,7);
plot(a4)
subplot(5,2,8);
plot(d4)
subplot(5,2,9);
plot(a5)
subplot(5,2,10);
plot(d5)figure;
shift_fa5 = abs(fft(a5,2000)); % 傅里叶变换
shift_fd5 = abs(fft(d5,2000)); % 傅里叶变换
shift_fd4 = abs(fft(d4,2000)); % 傅里叶变换
shift_fd3 = abs(fft(d3,2000)); % 傅里叶变换
shift_fd2 = abs(fft(d2,2000)); % 傅里叶变换
shift_fd1 = abs(fft(d1,2000)); % 傅里叶变换
plot(shift_fa5(2:1001),'r');title('shift fft transform');hold on;
plot(shift_fd5(2:1001),'g');title('shift fft transform');hold on;
plot(shift_fd4(2:1001),'b');title('shift fft transform');hold on;
plot(shift_fd3(2:1001),'y');title('shift fft transform');hold on;
plot(shift_fd2(2:1001),'r');title('shift fft transform');hold on;
plot(shift_fd1(2:1001),'g');title('shift fft transform');hold on;
小波分析是一种基于频率逐步二分的方法,信号x=a1+d1;a1、d1分量频率成分各为一半;a1=a2+d2,a2、d2分量频率成分各为一半,以此类推。
详细了解请科学上网
小波分解得到各成分的频率分布为:
5. 时域特征值
% 均方根、峰值因子、脉冲因子、裕度因子、峭度因子、波形因子和偏度
clear; clc;
% x=0:0.1:2*pi;
% y=sin(x); %信号
y = normrnd(0,2,2000,1); % 生成平均值为0,标准差为2的1000*1的矩阵
figure;
plot(y);
ma = max(y); % 最大值
mi = min(y); % 最小值
me = mean(y); % 平均值
pk = ma-mi; % 峰-峰值
av = mean(abs(y)); % 绝对值的平均值(整流平均值)
va = var(y); % 方差
st = std(y); % 标准差
ku = kurtosis(y); % 峭度
sk = skewness(y); % 偏度
rm = rms(y); % 均方根
S = rm/av; % 波形因子
C = pk/rm; % 峰值因子
I = pk/av; % 脉冲因子
xr = mean(sqrt(abs(y)))^2;
L = pk/xr; % 裕度因子
信号频谱分析、功率谱分析、倒谱分析、小波分析相关推荐
- 语音信号的同态处理、倒谱分析和Mel频率倒谱系数
1 同态处理 信号的同态处理也称同态滤波.大概步骤为: f(x,y)→ln→DFT→H(u,v)→(DFT)-1→exp→g(x,y) 虽然,一般用于图像处理.但是,博主将同态滤波用于语音信号的滤波. ...
- 语音信号同态与复倒谱分析
目录 1. 同态分析 2. 复倒谱与倒谱分析 3. 同态声码器 1. 同态分析 同态信号处理也称为同态滤波,其目的是将卷积关系变换为求和关系的分离处理.同态滤波是非线性滤波但服从广义叠加原理,即几种不 ...
- matlab频谱分析中振幅的物理意义,对速度信号进行傅里叶谱分析之后,其纵坐标对应的幅值的物理意义是什么?是速度,还是振幅...
横坐标是频率,2113纵坐标是对应5261频率成分的幅度.对4102速度信号进行傅里叶谱分析之后,纵1653坐标表示的是不同加速度的幅度.傅立叶原理表明:任何连续测量的时序或信号,都可以表示为不同频率 ...
- 连续语音信号的短时倒谱分析及其参数用途
文章目录 前言 基本概念 倒谱和倒谱参数 1.倒谱 2.复倒谱 3.倒谱计算 (1).由声门激励信号提取基音周期 (2).由声道冲激响应估算共振峰 4.倒谱的频谱 5.倒谱距离 6.Mel频率倒谱系数 ...
- matlab实验——信号和噪声产生及其功率谱分析
文章目录 前言 一.高斯白噪声自相关函数及功率谱 1.matlab代码 2.运行结果 二.均匀白噪声自相关函数及功率谱 1.matlab代码 2.运行结果 三.正弦波与高斯白噪声叠加 1.matlab ...
- 信号归一化功率_UE低发射功率余量分析
1.功率余量基础原理 1.1 功率余量PH 云南 (图1) 1.2 功率余量报告PHR PHR,全称是Power Headroom Report,中文为功率余量报告,即UE向网侧报告功率余量的过程.这 ...
- 倒谱分析 matlab,倒谱分析
看到一个倒谱分析的方法,网上查看到一些资料,介绍的挺不错的,分享一下: 这个是东南大学机械工程系的网上课件,略显粗糙,可以将就看 倒谱分析 用倒频谱诊断齿轮故障 分离信息通道对信号的影响 这个是一个新 ...
- 功率谱分析笔记-------脑电相关
1:功率谱分析的方法介绍 功率谱分析的方法大致可以分为两大类:第一类是经典的功率谱计算方法,第二类是现代功率谱计算方法,如图1所示.其中第一类经典功率谱分析方法,又可以分为直接法.间接法和改进的直接法 ...
- 语音信号的梅尔频率倒谱系数(MFCC)的原理讲解及python实现
梅尔倒谱系数(MFCC) 梅尔倒谱系数(Mel-scale FrequencyCepstral Coefficients,简称MFCC).依据人的听觉实验结果来分析语音的频谱, MFCC分析依据的听觉 ...
最新文章
- jstat 内存泄漏_一次Java内存泄漏的排查!要了自己的老命!
- php魔术方法例子,PHP中魔术方法的使用举例
- linux lamp实验报告,我的LAMP过程
- Now that half a year has passed
- 如何安装Pycharm和汉化包(包括安装软件,无广告)
- 史上最全的vue.js源码解析(一)
- 洛谷P3373线段树2
- 详解input value属性
- ML-Agents案例之双人足球
- WiFi手机可以连接,电脑上也能连接,可以微信聊天,但是不能浏览器上网怎么解决?
- CNVD原创漏洞证书总结
- 图像处理+帧差法实现对车辆的识别框选
- css3 太极动画,纯css实现太极阴阳鱼动画
- 看不懂没关系, 知道厉害就行了! 中科大俩教授11年解了两道数学难题
- 【Apache】Apache服务器配置
- 2020软考 信息安全工程师(第二版)学习总结【一】
- GEO数据库架构介绍
- 【书单】推荐书籍搜集
- python list遍历删除,对python list 遍历删除的正确方法详解
- linux cmds - ldconfig
热门文章
- [趣味分享]纯手工制作的滑翔机
- Selenium自动化测试的【投资回报率】还能这样计算
- 关于自定义tabbar项的问题
- GRS认证产销链监管.TC交易证书的常见问题
- Debug [Error] too many initializers for ‘<anonymous struct>‘
- HackThisSite(Chanllenages - Basic)
- m6c2g核心板使用笔记
- 抓取网页并保存静态资源
- Facebook,Twitter和实时聊天按钮:创建精通社交的WordPress网站
- coco数据集大小分类_COCO 数据集使用说明书