信号的频谱分析希尔伯特变换分析matlab仿真,含仿真操作录像

% function []=demo7() function []=InstantaneousCharactorUsingHilbertAndWavLet() %函数功能： % Hilbert方法计算信号的瞬时特征 Vs 计算小波系数方法提取信号的瞬时特征 %实验证明： % 1、小波方法提取MFSK信号的瞬时频率特征效果优于Hilbert方法计算信号的瞬时频率特征 % 2、但是，对于实际信号可能含噪较大或不是gaussian噪声，致使此法效果不理想；因此， % 应用小波处理前需进行降噪处理；然而，对于MFSK信号可能含有高频成分，用小波去噪 % 可能会使原信号发生不可逆恢复的破坏。 %码元速率为：瞬时频率峰值间隔的倒数或延迟自相关函数谱的峰值间隔 %Copyright@2009-2011 by hank(SiChuan University) %% close all;clear all;clc; %% %构造原始的4fsk信号 fs=1e4;dt=1/fs;bit=40;N=300;%波特率小于载频 % [Sig,tIndx]=fsk4Signal(dt,N,bit);%4fsk [Sig,tIndx]=fsk2Signal(dt,N,bit);%2fsk % [Sig,tIndx]=bpskSignal(dt,N,bit);%bpsk % [Sig,tIndx]=qpskSignal(dt,N,bit);%qpsk % [Sig,tIndx]=epskSignal(dt,N,bit);%8psk figure;subplot(221);plot(tIndx,Sig);title('纯净信号');axis([0 tIndx(end) -2 2]);%原始信号 SigR=awgn(Sig,20); subplot(222);plot(tIndx,SigR);title('染噪信号');axis([0 tIndx(end) -2 2]);%原始信号 nfft=4096;f=0:fs/nfft:fs/2-fs/nfft;index=1:nfft/2;%fft点数需要比载频大 SigRfft=abs(fft(SigR,nfft));Sigfft=abs(fft(Sig,nfft)); subplot(223);plot(f,Sigfft(index));title('纯净信号频谱'); subplot(224);plot(f,SigRfft(index));title('染噪信号频谱'); OddityDectation(SigR,tIndx);%检测信号的突变点 %% tau=0; nsamp=nfft;overlap=60;Window=hamming(nsamp); [P1,P2]=WelchPow(Sig',tau,nsamp,overlap,Window,nfft); figure;subplot(211);plot(f,10*log10(P1(index)));title('延迟功率谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); subplot(212);plot(f,10*log10(P2(index)));title('延迟自相关倒谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); %% %将原始信号分两路，一路进行hilbert变换后平方，另一路直接平方 %将两路信号叠加后求自相关，再进行fft进行码元速率估计(一次功率谱的谱峰间隔即为码速率) %当snr=0db时尚准(仅对psk信号有效) tau=0;Sigh=hilbert(Sig);SigUp=Sigh.*Sigh; SigDn=Sig.*Sig;SigMix=SigUp+SigDn; nsamp=nfft;overlap=60;Window=hamming(nsamp); [P1,P2]=WelchPow(SigMix',tau,nsamp,overlap,Window,nfft); figure;subplot(211);plot(f,10*log10(P1(index)));title('延迟功率谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); subplot(212);plot(f,10*log10(P2(index)));title('延迟自相关倒谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); %% %hilbert变换 SigH=hilbert(SigR);SigIm=imag(SigH); % figure;subplot(211);plot(tIndx,real(SigH));title('hilbert变换后信号的实部'); % subplot(212);plot(tIndx,imag(SigH));title('hilbert变换后信号的虚部'); %% %计算瞬时特征 [Amp,Ang,Freq]=IQinstantHeft(SigR,SigIm,fs); figure;subplot(311);plot(tIndx,Amp);title('瞬时幅度(hilbert法)'); subplot(312);plot(tIndx,Ang);title('瞬时相位(hilbert法)'); subplot(313);plot(tIndx,Freq);title('瞬时频率(hilbert法)'); %% nsamp=nfft;overlap=60;Window=hamming(nsamp); [P1,P2]=WelchPow(Amp',tau,nsamp,overlap,Window,nfft); figure;subplot(211);plot(f,10*log10(P1(index)));title('延迟功率谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); subplot(212);plot(f,10*log10(P2(index)));title('延迟自相关倒谱(Welch法)'); xlabel('频率/Hz');ylabel('幅度/db'); %% %提取小波系数 [AmpCoef1,AmpCoef2]=WaveLetCoefsExtra(SigR); figure;subplot(211);plot(tIndx,AmpCoef1);title('小波系数(功率归一化)'); subplot(212);plot(tIndx,AmpCoef2);title('小波系数(幅度归一化)'); %% %小波降噪后进行重复上述处理步骤 [SigDns1,SigDns2]=WaveLetNoiseRemove(AmpCoef1'); figure;subplot(211);plot(SigDns1);title('去噪后的小波系数(功率归一化)'); subplot(212);plot(SigDns2);title('去噪后的小波系数(功率归一化)'); [SigDns3,SigDns4]=WaveLetNoiseRemove(AmpCoef2'); figure;subplot(211);plot(SigDns3);title('去噪后的小波系数(幅度归一化)'); subplot(212);plot(SigDns4);title('去噪后的小波系数(幅度归一化)'); %% %降噪后小波系数的0延迟平方函数谱 SigDns3fft=abs(fft(SigDns4.*SigDns4,nfft)); figure;plot(f,10*log10((SigDns3fft(index)))); % SigDnsH=hilbert(SigDns1); % SigR2=hilbert(SigDnsH);SigIm2=imag(SigDnsH); % figure;subplot(211);plot(real(SigDnsH));title('hilbert变换后信号的实部'); % subplot(212);plot(imag(SigDnsH));title('hilbert变换后信号的虚部'); % [Amp1,Ang1,Freq1]=IQinstantHeft(SigR2,SigIm2,fs); % figure;subplot(311);plot(Amp1); % subplot(312);plot(Ang1);subplot(313);plot(Freq1); %% function [Sig,tIndx]=fsk2Signal(dt,N,bit) %函数功能：仿真一个2fsk信号 %输入参数： % dt:采样间隔 N:码元个数 bit:波特率 %输出参数： % Sig:产生的信号 tIndx：时间尺度 %% tb=1/bit;%码元周期 t0=dt:dt:tb; f=[1e3 1.8e3 ];%载频必须大于码速率 t=zeros(N,length(t0));y=zeros(N,length(t0)); for k=1:N t(k,:)=(k-1)*tb+t0; flag=mod(k,2); y(k,:)=cos(2*pi*f(flag+1)*t(1,:)); end [m,n]=size(y); Sig=reshape(y',1,m*n); tIndx=reshape(t',1,m*n); %% function [Sig,tIndx]=fsk4Signal(dt,N,bit) %函数功能：仿真一个4fsk信号 %输入参数： % dt:采样间隔 N:码元个数 bit:波特率 %输出参数： % Sig:产生的信号 tIndx：时间尺度 %% tb=1/bit;%码元周期 t0=dt:dt:tb; f=[1e3 1.8e3 2.5e3 3.3e3];%载频必须大于码速率 t=zeros(N,length(t0));y=zeros(N,length(t0)); for k=1:N t(k,:)=(k-1)*tb+t0; flag=mod(k,4); y(k,:)=cos(2*pi*f(flag+1)*t(1,:)); end [m,n]=size(y); Sig=reshape(y',1,m*n); tIndx=reshape(t',1,m*n); %% function [Sig,tIndx]=bpskSignal(dt,N,bit) %函数功能：仿真一个bpsk信号 %输入参数： % dt:采样间隔 N:码元个数 bit:波特率 %输出参数： % Sig:产生的信号 tIndx：时间尺度 %% tb=1/bit;%码元周期 t0=dt:dt:tb; f=1e3 ;%载频必须大于码速率 t=zeros(N,length(t0));y=zeros(N,length(t0));fai=[0 pi]; for k=1:N t(k,:)=(k-1)*tb+t0; flag=mod(k,2); y(k,:)=cos(2*pi*f*t(1,:)+fai(flag+1)); end [m,n]=size(y); Sig=reshape(y',1,m*n); tIndx=reshape(t',1,m*n); %% function [Sig,tIndx]=qpskSignal(dt,N,bit) %函数功能：仿真一个qpsk信号 %输入参数： % dt:采样间隔 N:码元个数 bit:波特率 %输出参数： % Sig:产生的信号 tIndx：时间尺度 %% tb=1/bit;%码元周期 t0=dt:dt:tb; f=1e3 ;%载频必须大于码速率 t=zeros(N,length(t0));y=zeros(N,length(t0));fai=[pi/4 3*pi/4 5*pi/4 7*pi/4]; for k=1:N t(k,:)=(k-1)*tb+t0; flag=mod(k,4); y(k,:)=cos(2*pi*f*t(1,:)+fai(flag+1)); end [m,n]=size(y); Sig=reshape(y',1,m*n); tIndx=reshape(t',1,m*n); %% function [Sig,tIndx]=epskSignal(dt,N,bit) %函数功能：仿真一个8psk信号 %输入参数： % dt:采样间隔 N:码元个数 bit:波特率 %输出参数： % Sig:产生的信号 tIndx：时间尺度 %% tb=1/bit;%码元周期 t0=dt:dt:tb; f=1e3 ;%载频必须大于码速率 t=zeros(N,length(t0));y=zeros(N,length(t0)); fai=[pi/8 3*pi/8 5*pi/8 7*pi/8 9*pi/8 11*pi/8 13*pi/8 15*pi/8]; for k=1:N t(k,:)=(k-1)*tb+t0; flag=mod(k,8); y(k,:)=cos(2*pi*f*t(1,:)+fai(flag+1)); end [m,n]=size(y); Sig=reshape(y',1,m*n); tIndx=reshape(t',1,m*n); %% function [SigDns1,SigDns2]=WaveLetNoiseRemove(Sig) %函数功能：小波变换对信号进行降噪 %输入参数： % Sig:输入信号 %输出参数： % SigDns：降噪后的信号 %% % [c,l]=wavedec(Sig',2,'db5'); % sigma=wnoisest(c,l,1); % alpha=2; % thr=wbmpen(c,l,sigma,alpha); % keepapp=1; % SigDns=wdencmp('gbl',c,l,'db5',2,thr,'s',keepapp); %% %根据信号计算噪声强度，给出全局阈值 [thr,sorh,keepapp]=ddencmp('den','wv',Sig); %根据全局阈值对信号降噪 SigDns1=wdencmp('gbl',Sig,'sym4',4,thr,sorh,keepapp); %用sym4小波对信号做4层分解 [c,l]=wavedec(Sig,4,'sym4'); %得到每层的分层阈值 [thr1,nkeep]=wdcbm(c,l,2); %根据分层阈值使用软阈值方法对信号进行降噪 [SigDns2,cxd,lxd,perf0,perf12]=wdencmp('lvd',c,l,'haar',4,thr1,'s'); %% function [Amp,Ang,Freq]=IQinstantHeft(I,Q,fs) %函数功能： % 分析IQ分量信号的瞬时特征信息 %输入参数: % I/Q：同相分量/正交分量 fs：采样频率 %输出参数: % Amp: 瞬时幅度 % Ang: 真正的瞬时相位 % Freq: 瞬时频率 %% N=length(I); H=I+j*Q;Amp=abs(H);Ang=angle(H); AngAmen=zeros(N,1); Freq=zeros(N-1,1); for k=2:N Freq(k)=1/(2*pi)*fs*( ( I(k-1)*Q(k) )- ( I(k)*Q(k-1) ) )/( I(k)^2+Q(k)^2); end %% function [AmpCoef1,AmpCoef2]=WaveLetCoefsExtra(Sig) %函数功能：小波变换提取小波系数 %输入参数： % Sig:输入信号 %