(1)傅里叶变换和希尔伯特变换及相关函数在信号处理中的MATLAB仿真

     这次博文主要内容是对离散傅里叶变换和希尔伯特变换进行简单介绍。在MATLAB中仿真出一些例子，让我们学习时对理论知识有一个形象的理性认识，这样理解起来比较深刻。

    内容概要：利用离散傅里叶变换将加噪的调制信号变换到频域，用去除高频的高斯白噪声干扰的方法进行降噪。然后利用希尔伯特变换求得调制信号的解析信号，根据解析式得到调制信号的瞬时参数：瞬时幅度，瞬时频率，瞬时相位。最后用自相关函数检测调制信号的码元速率。

 

一．理论分析（理论部分这里博主略去，很多公式不方便录入，学习者自己找，很多）

1.离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。

2.希尔伯特变换

希尔伯特变换与傅里叶变换不同，它不是把信号从时间域变换到另外的域，而是把信号从时域仍然变换到时域。

3.自相关函数

二．针对问题和方法

1.信号降噪：现实中的信号一般都受到噪声的干扰，而导致在提取信号特征信息的造成误差，所以在处理信号之前先进行降噪。一个含噪声的一维信号模型可表示为如下形式：

S(k)=f(k)+e(k)

S(k)为含噪信号，f(k)是有用信号，e(k)是噪声信号。这里假定噪声是高斯白噪声，频谱一般分布为整个频域。

 

20db信噪比4psk的频谱

实际工程中f(k)通常为低频信号或者频谱范围分布有限的信号。因此，通过离散傅里叶变换得到信号的频谱以后，可以把有用的频谱之外的噪声频谱去除，然后再经过离散傅里叶逆变换就可以得到降噪的信号。

 

20db信噪比4psk去除一段噪声的频谱

 

2.提取调制信号瞬时参数：利用信号的解析函数的特性计算调制信号的特征。

用原信号和其希尔伯特变换表示的复函数称作这个原实信号的解析函数

z(tn)=u(n)+jv(n)

采用频域的方法计算采样信号u(n)的希尔伯特变换v(n)。

 

 

3.利用自相关函数求调制信号码元速率

         调制信号在相位突变的时候，由解析函数求得的瞬时频率在码元跳变处会频率突变。利用这一点来检测调制信号的码元速率。先把瞬时频率信号减去瞬时频率的平均值，使得频率中心尽量归零，接着对处理后的瞬时频率信号进行自相关，自相关后的信号每隔一段采样点就会出现突变，这些突变的最小间隔时间就是码元速率Rs。

 

三.MATLAB仿真

调制信号的瞬时参数在很多方面都有重要应用，下面针对数字相位调制信号4psk（加高斯白噪声），首先进行降噪，然后根据解析信号的求得瞬时参数。仿真如下：

         仿真信号s(n)采用4psk离散调制信号，载波选择fs=1kHz，码元速率fd=100Hz，采样频率选择fs=10kHz，叠加信噪比20db的高斯白噪声。

 

信噪比为20db的4psk信号（未降噪）瞬时参数提取结果：

 

 

信噪比为20db的4psk信号（降噪后）瞬时参数提取结果：

 

由这两幅仿真图看出，信号降噪前后有明显的效果。由解析信号提取的瞬时幅度在1附近变化，瞬时频率在1kHz附近变化，在相位突变处瞬时频率发生剧烈变化。调制信号相位对应的码元xi=[0 1 2 3 3 2 1 0 2 3 0 1]

 

瞬时频率信号归零后的自相关的绝对值：

   由图可得，基本每隔Nt=100个采样点就会出现一个自相关函数值的剧烈增大，而信号的采样频率fs=10000Hz，所以码元速率Rs=Nt/fs=0.01b/s.调制的码元频率fd=100Hz，理论上码元速率是1/fd,即0.01b/s。码元速率的检测值和理论值接近。

 

主m程序

clear all

 

M=4;

fc=1000;

fd=100;

fs=10000;

xi=[0 1 2 3 3 2 1 0 2 3 0 1];

[u,t]=dmod(xi,fc,fd,fs,'psk',M);

sn0=awgn(u,20,'measured','db');%?ó??éù

sn=fftlvbo(sn0);

z=hilbert(sn);

fai1=angle(z);

p=FXXXWTQ3(z,fc,fs);

 

fi=f(z,fc,fs);

z=abs(z);

figure(1)

subplot(2,2,1);

plot(t,sn);

title('?ó??μ÷??D?o?sn=f(t)')

 

subplot(2,2,2);

plot(t,z);

title('μ÷??D?o??2ê±·ù?è|z|')

subplot(2,2,3);

plot(t,fi);

title('μ÷??D?o??2ê±?μ?êfi')

 

subplot(2,2,4);

dt=1/fs;

Np=length(p);

tp=[1:Np]*dt;

plot(tp,p);

title('μ÷??D?o??2ê±?à??p￡¨?à??è￥??μt￡?')

 

fi1=fi-mean(fi);

r = xcorr(fi1);

figure(3)

plot(abs(r))

 

降噪m程序

function u=fftlvbo(sn)

f=fft(sn);

L=length(f);

for i=270:930

    if(abs(f(i))<20)

        f(i)=0;

    end

end

u=ifft(f);

figure(2)

plot(abs(f))

 

解析信号构建m文件

U=fft(u,N);%à?ó??￡??2?ì?±???￡?′ó?μóò±???μ?μ?v￡¨t￡?

k=1:N/2;%V(f)=-j*sgn(f)*U(f)

V(k)=-j*U(k);

k=(N+2)/2:N;

V(k)=j*U(k);

v=ifft(V,N);%?μà?ò???±???μ?μ?v￡¨t￡?

z=u+j*v;%?a??D?o?

作者；西电大宝