查看: 1850|回复: 16|关注: 0
关于用HHT求信号的瞬时频率和时频谱
关注者： 2
设信号是y=5xin(2*pi*10t)+5*sin(2*pi*35t)，采样频率为2048Hz
%用HHT求信号的是频谱和边际谱
function HHT
%fft默认计算的信号时从0开始的
t=linspace(1,2,N);
deta=t(2)-t(1);
x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
%计算每个IMF分量及最后一个剩余分量residual与原信号的相关性
[m,n]=size(c);
for i=1:m;
& & a=corrcoef(c(i,:),z);%corrcoef函数是计算两组数据的相关系数，corrcoef(x,y)表示x序列和y序列的相关系数，得到的是2*2矩阵，对角元素为自相关系数，非对角元素为互相关系数
& & xg(i)=a(1,2);
for i=1:m-1;
& & %计算各IMF的方差贡献率
& & %定义：方差=平方的均值-均值的平方，即D(x)=E(x^2)-[E(x)]^2
& & %平方的均值 imf2p=mean(c(i,:).^2,2)
& & %均值的平方 imfp2=mean(c(i,:),2).^2
%各IMF的方差
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
mmse=sum(mse);%mmse为所有的方差和
for i=1:m-1;
& & mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
& & %方差百分比，也就是方差贡献率
& & mseb(i)=mse(i)/mmse*100;
%画出每个IMF分量及最后一个剩余分量residual的图形
for i=1:m-1;
& & disp(['imf',int2str(i)]);%int2str(i)将数字i四舍五入转变为字符，若i=3.3，则str1=imf(3);
& & disp([mse(i),mseb(i)]);%显示各个IMF的方差和其方差贡献率
subplot(m+1,1,1);%有m个IMF和1个剩余量，故划分m+1行和1列
%画原始信号的图形
plot(t,z);
set(gca,'fontname','Times New Roman');%设定gca的字体为times new roman
set(gca,'fontsize',14.0);%设定gca的字号为14.0
ylabel(['signal','amplitude']);%y轴设定为原信号的幅值
%从第二行开始画imf1~imfm的图形
for i=1:m-1;
& & subplot(m+1,1,i+1);%画地IMF(I+1)的图形
& & set(gcf,'color','w');%设定gca的颜色为白色
& & plot(t,c(i,:),'k');%画时间t-imf(i)的图形
& & set(gca,'fontname','times new roman');
& & set(gca,'fontsize',14.0);
& & ylabel(['imf',int2str(i)]);%y轴为imf
%画剩余量的图（从第m+1行开始）
subplot(m+1,1,m+1);
set(gcf,'color','w');
plot(t,c(m,:),'w');
set(gca,'fontname','times new roman');
set(gca,'fontsize',14.0);
ylabel(['r',int2str(m-1)]);%y轴为剩余量r
%画出每个IMF分量及剩余量residual的幅频曲线
subplot(m+1,1,1)
set(gcf,'color','w');
%[f,z]=fftfenxi(t,z);%对时间-幅值矩阵做fft求的原始信号的频谱图
plot(f,z,'k');
set(gca,'fontname','times new roman');
set(gca,'fontsize',14.0);
ylabel(['initial signal',int2str(m-1),'amplitude']);%y轴为原始信号的幅值
%画各个IMF(i)的频谱图
for i=1:m-1;
& & subplot(m+1,1,i+1);
& & set(gcf,'color','w');
& & [f,z]=fftfenxi(t,c(i,:));
& & plot(f,z,'k');
& & set(gca,'fontname','times new roman');
& & set(gca,'fontsize',14.0);
& & ylabel(['imf',int2str(i),'amplitude']);%y轴为imf(i)的幅值
%画剩余量的频谱图
subplot(m+1,1,m+1);
set(gcf,'color','w');
[f,z]=fftfenxi(t,c(m,:));
plot(f,z,'k');
set(gca,'fontname','times new roman');
set(gca,'fontsize','14.0');
ylabel(['r',int2str(m-1),'amplitude']);
hx=hilbert(z);
xr=real(hx);
xi=imag(hx);
%计算瞬时政府sz
sz=sqrt(xr.^2+xi.^2);
%计算瞬时相位sx
sx=angle(hx);
%计算瞬时频率sp
dt=diff(t);
dx=diff(sx);
plot(t(1:N-1),sp);
title('瞬时频率');
提示出现的错误是：
??? Error using ==& HHT&fftfenxi
Too many input arguments.
Error in ==& HHT at 85
& & [f,z]=fftfenxi(t,c(i,:));
我是新手，很多都还不懂，有哪位高手能否解答一下呀，十分感谢。。
关注者： 7
函数的输入变量多了。
关注者： 2
函数的输入变量多了。
我输入的是一序列，由不同的时间，采样的点。那这个要怎么改呢？
关注者： 7
我输入的是一序列，由不同的时间，采样的点。那这个要怎么改呢？
LZ最好看看你用的函数变量说明。
关注者： 2
LZ最好看看你用的函数变量说明。%用HHT求信号的是频谱和边际谱
function HHT
N=2048;
%fft默认计算的信号时从0开始的
t=linspace(1,2,N);
deta=t(2)-t(1);
fs=1/
x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);
z=x;
c=emd(z);
%计算每个IMF分量及最后一个剩余分量residual与原信号的相关性
[m,n]=size(c);
for i=1:m;
& & a=corrcoef(c(i,:),z);%corrcoef函数是计算两组数据的相关系数，corrcoef(x,y)表示x序列和y序列的相关系数，得到的是2*2矩阵，对角元素为自相关系数，非对角元素为互相关系数
& & xg(i)=a(1,2);
end
for i=1:4;
& & %计算各IMF的方差贡献率
& & %定义：方差=平方的均值-均值的平方，即D(x)=E(x^2)-[E(x)]^2
& & %平方的均值 imf2p=mean(c(i,:).^2,2)
& & %均值的平方 imfp2=mean(c(i,:),2).^2
%各IMF的方差
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
end
mmse=sum(mse);%mmse为所有的方差和
for i=1:4;
& & mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
& & %方差百分比，也就是方差贡献率
& & mseb(i)=mse(i)/mmse*100;
end
%画出每个IMF分量及最后一个剩余分量residual的图形
figure(1)
for i=1:4;
& & disp(['imf',int2str(i)]);%int2str(i)将数字i四舍五入转变为字符，若i=3.3，则str1=imf(3);
& & disp([mse(i),mseb(i)]);%显示各个IMF的方差和其方差贡献率
end
subplot(m+1,1,1);%有m个IMF和1个剩余量，故划分m+1行和1列
%画原始信号的图形
plot(t,z);
set(gca,'fontname','Times New Roman');%设定gca的字体为times new roman
set(gca,'fontsize',14.0);%设定gca的字号为14.0
ylabel(['signal','amplitude']);%y轴设定为原信号的幅值
%从第二行开始画imf1~imfm的图形
for i=1:4;
& & subplot(6,1,i+1);%画地IMF(I+1)的图形
& & set(gcf,'color','w');%设定gca的颜色为白色
& & plot(t,c(i,:),'k');%画时间t-imf(i)的图形
& & set(gca,'fontname','times new roman');
& & set(gca,'fontsize',14.0);
& & ylabel(['imf',int2str(i)]);%y轴为imf
end
%画剩余量的图（从第m+1行开始）
subplot(6,1,6);
set(gcf,'color','w');
plot(t,c(3,:),'w');
set(gca,'fontname','times new roman');
set(gca,'fontsize',14.0);
ylabel(['r',int2str(3-1)]);%y轴为剩余量r
%画出每个IMF分量及剩余量residual的幅频曲线
f=m/N*
figure(2)
subplot(6,1,1)
set(gcf,'color','w');
%[f,z]=fftfenxi(t,z);%对时间-幅值矩阵做fft求的原始信号的频谱图
plot(f,z,'k');
set(gca,'fontname','times new roman');
set(gca,'fontsize',14.0);
ylabel(['initial signal',int2str(3-1),'amplitude']);%y轴为原始信号的幅值
%画各个IMF(i)的频谱图
for i=1:4;
& & subplot(6,1,i+1);
& & set(gcf,'color','w');
& & [f,z]=fftfenxi(t,c(i,:));
& & plot(f,z,'k');
& & set(gca,'fontname','times new roman');
& & set(gca,'fontsize',14.0);
& & ylabel(['imf',int2str(i),'amplitude']);%y轴为imf(i)的幅值
end
%画剩余量的频谱图
subplot(6,1,6);
set(gcf,'color','w');
[f,z]=fftfenxi(t,c(3,:));
plot(f,z,'k');
set(gca,'fontname','times new roman');
set(gca,'fontsize',14.0);
ylabel(['r',int2str(3-1),'amplitude']);
hx=hilbert(z);
xr=real(hx);
xi=imag(hx);
%计算瞬时政府sz
sz=sqrt(xr.^2+xi.^2);
%计算瞬时相位sx
sx=angle(hx);
%计算瞬时频率sp
dt=diff(t);
dx=diff(sx);
sp=dx./
figure(6)
plot(t(1:N-1),sp);
title('瞬时频率');
end
function hhspectrum
[A,fa,tt]=hhspectrum(c);%A为IMF的幅值，fa为瞬时频率，tt为截断时间瞬间
[E,tt1]=toimage(A,fa,tt,length(tt));%E为2维的频谱图，tt1为图像的瞬时时间，length(tt)为信号的长度
figure(3)
disp_hhs(E,tt1);%用二维图显示HHT时频谱，E是求得的HHT谱
pause
figure(4)
for i=1:size(c,1);
& & faa=fa(i,:);%faa为由各个1~n个IMF对应的瞬时频率对应的瞬时频率矩阵
& & [FA.TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图
& & surf(FA,TT1,E)%surf是产生三维曲面（色）图，surf(x,Y,Z)能够产生有X、Y、Z指定的有色参数化曲面，即三维有色图
& & title('HHT时频谱三维显示');
& & hold on
end
hold off
E=flipud(E);%flipud是实现矩阵的上下翻转，而fliplr是实现矩阵的左右翻转
for k=1:siza(E,1);
& & bjp(k)=sum(E(k,:))*1/%边际谱
end
f=(1:N-2)/N*(fs/2); %频率序列
figure(5)
plot(f,bjp);
xlabel('频率/Hz');
ylabel('信号幅值');
title('信号边际谱');%要求边际谱必须先对信号进行EMD分解
end
function[A,f,tt]=hhspectrum(x,t,l,aff)
error(nargin(1,4,nargin));%nargin函数参数个数
if nargin&2
& & t=1:size(x,2);
end
if nargin&3
& & l=1;
end
if nargin&4
& & aff=0;
end
if min(size(x))==1
& & if size(x,2)==1;
& && &&&x=x';
& && &&&if nargin&2
& && && && &t=1:size(x,2);
& && &&&end
& & end
& & Nmodes=1;
else
& & Nmodes=size(x,1);
end
lt=length(t);
tt=t((l+1):(lt_1));
for i=1:Nmodes
& & an(i,:)=hilbert(x(i,:)')';
& & f(i,:)=instfreq(an(i,:)',tt,l)';
& & A=abs(an(:,l+1:end-1));
& & if aff
& && &&&disprog(i,Nmodes,max(Nmodes,100));
& & end
end
end
function&&disp_hhs(im,t,inf)
%disp_hhs(im,t,inf)
%displays in a new figure the spectrum contained in log(用log表示振幅）
%amplitudes in log(用Log表示振幅）
%输入：
%im:图像矩阵（如toimage的输出）
%t:瞬时时间（如toimage输出）
%inf:用dB(分贝）的动态范围，其默认值：inf=-20
figure
colormap(bone)%将颜色从黑变白
colormap(1-colormap);
if nargin==-1;
& & inf=-20;
& & t=1:size(in,2);
end
if nargin==2
& & if length(t)==1;
& && &&&inf=t;
& && &&&t=1:size(im,2);
& & else
& && &&&inf=-20;
& & end
end
if inf&=0
& & error('inf do it etre&0');
end
M=max(max(im));
im=log19(im/M+1e-300);
inf=inf/10;
imagesc(t,fliplr((1:size(im,1))/(2*size(im,1))),im,[inf,0]);%imagesc函数显示灰度图像
set(gca,'YDir',normal);
xlabel('时间');
ylabel('标准化频率');
tltle('希尔伯特黄谱Hilbert-Huang spectrum');
end
function fftfenxi
& && &&&[f,z]=fftfenxi(t,y);
& && &&&L=length(t);
& && &&&N=2.^nextpow2(L);
%fft默认计算的信号是从0开始的
t=linspace(t(1),t(L),N);
deta=t(2)-t(1);
m=0:N-1;
f=1./(N*deta)*m;%频率序列
%计算的Y就是x(t)的傅里叶变换值
%Y=exp(i*4*pi*f).*fft(y);
%将计算的频谱乘以exp（i*4*pi*f)得到频移后【-2,2】之间的频谱值
Y=fft(y);
z=sqrt(Y*conj(Y));
end复制代码我还是不是很懂，这是我的程序，你能否帮我看看。
本帖最后由 youyouma6262 于
14:49 编辑
缺少 FFTFENXI的程序，附件 中是程序，我也是从论坛里下载到的。
14:48 上传
点击文件名下载附件
738 Bytes, 下载次数: 69
关注者： 2
缺少 FFTFENXI的程序，附件 中是程序，我也是从论坛里下载到的。
我代码的最后一个函数有function fftfenxi函数呀？
你好，你这个问题解决了没有啊，我也遇到了同样的问题，不知道怎么改啊
关注者： 2
你好，你这个问题解决了没有啊，我也遇到了同样的问题，不知道怎么改啊
再加一个fftfenxin的函数就可以了。
再加一个fftfenxin的函数就可以了。
刚接触，不太懂哎，就比如下面的程序
function HHT
%fft默认计算的信号是从0开始的
t=linspace(1,2,N);deta=t(2)-t(1);fs=1/
x=sin(100*pi*t)+1.7*(heaviside(t-0.028)-heaviside(t-0.03))
%计算每个IMF分量及最后一个剩余分量residual与原始信号的相关性
[m,n]=size(c);
for i=1:m;
a=corrcoef(c(i,:),z);
xg(i)=a(1,2);
for i=1:m-1
%--------------------------------------------------------------------
%计算各IMF的方差贡献率
%定义：方差为平方的均值减去均值的平方
%均值的平方
%imfp2=mean(c(i,:),2).^2
%平方的均值
%imf2p=mean(c(i,:).^2,2)
%各个IMF的方
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
mmse=sum(mse);
for i=1:m-1
mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;
%方差百分比，也就是方差贡献率
mseb(i)=mse(i)/mmse*100;
%显示各个IMF的方差和贡献率
%画出每个IMF分量及最后一个剩余分量residual的图形
for i=1:m-1
disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
subplot(m+1,1,1)
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['signal','Amplitude'])
for i=1:m-1
subplot(m+1,1,i+1);
set(gcf,'color','w')
plot(t,c(i,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['imf',int2str(i)])
subplot(m+1,1,m+1);
set(gcf,'color','w')
plot(t,c(m,:),'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['r',int2str(m-1)])
%画出每个IMF分量及剩余分量residual的幅频曲线
subplot(m+1,1,1)
set(gcf,'color','w')
[f,z]=fftfenxi(t,z);
plot(f,z,'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['initial signal',int2str(m-1),'Amplitude'])
for i=1:m-1
subplot(m+1,1,i+1);
set(gcf,'color','w')
[f,z]=fftfenxi(t,c(i,:));
plot(f,z,'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['imf',int2str(i),'Amplitude'])
subplot(m+1,1,m+1);
set(gcf,'color','w')
[f,z]=fftfenxi(t,c(m,:));
plot(f,z,'k')
set(gca,'fontname','times New Roman')
set(gca,'fontsize',14.0)
ylabel(['r',int2str(m-1),'Amplitude'])
hx=hilbert(z);
xr=real(hx);xi=imag(hx);
%计算瞬时振幅
sz=sqrt(xr.^2+xi.^2);
%计算瞬时相位
sx=angle(hx);
%计算瞬时频率
dt=diff(t);
dx=diff(sx);
plot(t(1:N-1),sp)
title('瞬时频率')
应该怎么弄啊，谢谢啊，帮帮我！！
站长推荐 /1
Powered by把接收信号划分成为多个FFT窗口以减少互相关系数算术运算过程的路径寻找电路　 　
把接收信号划分成为多个FFT窗口以减少互相关系数算术运算过程的路径寻找电路　 　
申请专利号
专利申请日
把接收信号划分成为多个FFT窗口以减少互相关系数算术运算过程的路径寻找电路　 　
公开（公告）号
公开（公告）日
申请（专利权）
日本电气株式会社 　
日本东京都　
发明（设计）人
佐藤俊文 　
进入国家日期
专利代理机构
中国专利代理(香港)有限公司 　
一种路径寻找电路，接收的信号分为多个FFT窗口以减少对互相关系数算术运算处理。交错单元划分接收的信号为以一片时间间隙的rxd1,rxd2，和由重叠的FFT窗口选出该两个序列，和对由两个FFT单元选出的序列执行FFT。互功率谱计算单元确定FFT后的接收信号与存储在参考信号存储单元中的参考信号间的互功率谱。对于各FFT窗口由平均单元对互功率谱计算单元的输出进行平均，和由IFFT单元对平均的互功率谱执行IFFT。
权利要求书
1．用于检测路径定时的使用DS-CDMA通信方法的接收机的路径寻找电
路，所述路径定时是这样的定时，在发射端在这样的定时下扩展来自接收的无
线信号，所述电路包括：
用于对接收的无线电信号进行滤波和频率转换，以便将接收的无线电信号
转换成为基带信号的无线电接收单元；
以等于片速率N倍的取样速率取样基带信号、以便把基带信号转换成为
数字信号的A/D转换器；
互相关系数计算单元，包括把由所述A/D转换器数字化的基带信号重新
安排成为以片间隙取样的N个序列的交错装置，用来获取由交错装置用相互重
叠的预定时间长度的FFT窗口重新安排的N个接收信号序列和对选出的接收
信号序列执行快速傅利叶变换的N个快速傅利叶变换装置，用以存储通过选出
具有固定时间长度的FFT窗口的预定码序列和快速傅利叶变换选出的作为参考
信号的码序列而产生的信号序列的参考信号存储装置，用以为确定接收信号与
预定码序列间的互功率谱而为每个FFT窗口的确定由快速傅利叶变换装置快速
傅利叶变换的接收信号和存储在参考信号存储装置中的参考信号的共轭复数的
积的N个互功率谱计算装置，用以平均各个FFT窗口的互功率谱的N个互功
率谱平均装置，用以对由所述互功率谱平均装置平均的N个互功率谱进行反向
快速傅利叶变换以便把互功率谱转换成为N个互相关系数和输出该N个互相
关系数的N个反向快速傅利叶变换装置，和以产生和输出单个互相关系数的时
间顺序重新安排由所述各反向快速傅利叶变换装置输出的N个互相关系数的去
交错装置；
对由互相关系数计算单元以固定的时间周期输出的互相关系数进行平均的
互相关系数平均单元；和
从所述互相关系数平均单元平均的互相关系数中检测一个或多个峰值并将
得到峰值或每个峰值的定时作为路径定时输出的峰值检测单元。数字信号处理第四章9 线性卷积和线性相关的FFT算法_百度文库
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
文档贡献者
评价文档：
喜欢此文档的还喜欢
数字信号处理第四章9 线性卷积和线性相关的FFT算法
数​字​信​号​处​理
把文档贴到Blog、BBS或个人站等：
普通尺寸(450*500pix)
较大尺寸(630*500pix)
大小：327.50KB
登录百度文库，专享文档复制特权，财富值每天免费拿！
你可能喜欢信号相关系数 相关系数 相关系数公式 matlab 相关系数 excel 相关系数 pearson相关..
扫扫二维码，随身浏览文档
手机或平板扫扫即可继续访问
数字信号处理课程设计(两个序列之间的相关系数)
举报该文档为侵权文档。
举报该文档含有违规或不良信息。
反馈该文档无法正常浏览。
举报该文档为重复文档。
推荐理由：
将文档分享至：
分享完整地址
文档地址：
粘贴到BBS或博客
flash地址：
支持嵌入FLASH地址的网站使用
html代码：
&embed src='/DocinViewer-4.swf' width='100%' height='600' type=application/x-shockwave-flash ALLOWFULLSCREEN='true' ALLOWSCRIPTACCESS='always'&&/embed&
450px*300px480px*400px650px*490px
支持嵌入HTML代码的网站使用
您的内容已经提交成功
您所提交的内容需要审核后才能发布，请您等待！
3秒自动关闭窗口KARHUNENLOEVE变换及其在信号处理中的应用应用数学（APPLIED MATHEM..
扫扫二维码，随身浏览文档
手机或平板扫扫即可继续访问
KARHUNENLOEVE变换及其在信号处理中的应用
举报该文档为侵权文档。
举报该文档含有违规或不良信息。
反馈该文档无法正常浏览。
举报该文档为重复文档。
推荐理由：
将文档分享至：
分享完整地址
文档地址：
粘贴到BBS或博客
flash地址：
支持嵌入FLASH地址的网站使用
html代码：
&embed src='/DocinViewer-4.swf' width='100%' height='600' type=application/x-shockwave-flash ALLOWFULLSCREEN='true' ALLOWSCRIPTACCESS='always'&&/embed&
450px*300px480px*400px650px*490px
支持嵌入HTML代码的网站使用
您的内容已经提交成功
您所提交的内容需要审核后才能发布，请您等待！
3秒自动关闭窗口