用matlab处理信号,信号处理matlab程序示例

本人研究信号处理算法，利用matlab对信号做分析和处理。以下是我做的一个matlab程序实例，特贴出来共享，一起探讨学习。该程序绘制了仿真信号的FFT频谱和功率谱，并调用emd分析方法，绘制了emd分解结果，以及瞬时包络和瞬时频率图，还有Hilbert-Huang谱和Hilbert边际谱。

matlab程序代码如下：

clear all;

clc;

%给定参数

fs=1000;

N=1000;

n=0:N-1;

t=n/fs;

%构造仿真信号

x=2sin(2pi120t)+(1+sin(2pi5t)).sin(2pi50t+0.5cos(2pi5*t));

%emd分解

imf=emd(x);

%绘制时域波形

figure,plot(t,x);xlabel('时间

t/s');ylabel('幅值');title('仿真信号时域波形');

%绘制FFT频谱

nfft= 2^nextpow2(length(t));%找出大于y的个数的最大的2的指数值

ff=fs*(0:nfft/2-1)/nfft;�T变换后对应的频率的序列

fftx=fft(x,nfft);%求FFT变换

ps=fftx.*conj(fftx)/nfft;%conj()函数是求y函数的共轭复数，实数的共轭复数是他本身

figure;

subplot(211),plot(ff,abs(fftx(1:nfft/2))*2/nfft);ylabel('幅值');xlabel('频率');title('FFT频谱');

subplot(212),plot(ff,ps(1:nfft/2));ylabel('功率谱密度');xlabel('频率');title('信号功率谱');

%绘制emd分解结果

plot_imf(x,t,imf);

%求时频谱

[A,f,t]=hhspectrum(imf(1:end-1,

:));%对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号

%绘制瞬时包络图和瞬时频率图

figure;

subplot(221),plot(t/N,A(1,:));xlabel('时间

t/s');ylabel('幅值');title('imf1分量瞬时包络');

subplot(222),plot(t/N,f(1,:)*fs);xlabel('时间

t/s');ylabel('频率');title('imf1分量瞬时频率');

subplot(223),plot(t/N,A(2,:));xlabel('时间

t/s');ylabel('幅值');title('imf2分量瞬时包络');

subplot(224),plot(t/N,f(2,:)*fs);xlabel('时间

t/s');ylabel('频率');title('imf2分量瞬时频率');

%即时频图(用颜色表示第三维值的大小)和三维图(三维坐标系：时间，中心频率，振幅)

[E,t,Cenf]=toimage(A,f,t,length(t));%将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率 这里横轴为时间，纵轴为频率

%绘制Hilbert-Huang谱

figure;

set(gcf,'Color','w');

imagesc(t/N,[0,0.5*fs],E);

set(gca,'YDir','normal')

colormap('jet')

colorbar;

xlabel('时间 t/s');

ylabel('频率 f/Hz');axis([0 1 0 200])

title('Hilbert-Huang Spectrum') ;

%画出边际谱

%N=length(Cenf);%设置频率点数 %完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解

for k=1:size(E,1)

bjp(k)=sum(E(k,:))*1/fs;

end

figure;

plot(Cenf(1,:)*fs,bjp); %

作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率

xlabel('频率 f/Hz');

ylabel('幅值');

title('边际谱');

仿真结果如下：