阶跃函数的谱
邹俊杰
2016年11月14日
0 library
library(ggplot2)
library(stats)
library(base)
library(gridExtra)1 基础概念
傅里叶分析包含傅里叶级数与傅里叶变换。傅里叶级数用于对周期信号转换，傅里叶变换用于对非周期信号转换。但是对于不收敛信号，傅里叶变换无能为力，只能借助拉普拉斯变换。（主要用于计算微分方程）。而z变换则可以算作离散的拉普拉斯变换。（主要用于计算差分方程）。
1.1 傅里叶变换与拉普拉斯变换的关系
连续时间傅里叶变换的公式是： \[\int_{-\infty}^{\infty } x(t)e^{-jwt} dt,\] 这里的\(w\)是实数。
傅里叶变换要求时域信号绝对可积，即 \[\int_{-\infty }^{\infty } | x(t) | dt<\infty, \]
为了让不符合这个条件的信号，也能变换到频率域，我们给x(t)乘上一个指数函数\(e^{-\sigma t}\) ，\(\sigma\) 为任意实数。
可以发现，\(x(t)e^{-\sigma t}\) 这个函数，就满足了绝对可积的条件，即 \[\int_{-\infty }^{\infty } \left| x(t)e^{-\sigma t} \right| dt<\infty, \]
于是这个新函数的傅立叶变换就是： \[\int_{-\infty }^{\infty } x(t)e^{-\sigma t} e^{-j\omega t} dt,\]
化简得 \[\int_{-\infty }^{\infty } x(t)e^{-(\sigma +j\omega )t}\]
显然\(\sigma +j\omega\)是一个复数，我们把这个复数定义为一个新的变量——复频率，记为s。于是便得到了拉普拉斯变换的公式： \[\int_{-\infty }^{\infty } x(t)e^{-st} dt\]
拉普拉斯变换解决了不满足绝对可积条件的连续信号，变换到频率域的问题，同时也对“频率”的定义进行了扩充。
所以拉普拉斯变换与连续时间傅里叶变换的关系是：拉普拉斯变换将频率从实数推广为复数，因而傅里叶变换变成了拉普拉斯变换的一个特例。当s为纯虚数时，x(t)的拉普拉斯变换，即为x(t)的傅里叶变换。
另外，以上提到为双边拉氏变换，其收敛域可能为空。我们常用的是单边拉氏变换，他的收敛域必定非空。我们研究因果关系，主要关注 \(t>0\) 的系统响应，即便 \(t<0\) 网络内部有储能，但其影响可以用 \(t=0\) 时的初始条件表征。所以初始条件与 \(t\ge 0\) 的激励已经能够完全确定 \(t>0\) 的系统响应。
1.2 广义傅里叶变换
把原来不满足绝对可积的信号，通过乘上收敛因子取代极限的办法而求得的傅立叶变换，称为广义傅立叶变换。
1.3 功率谱与频谱，能量谱
功率谱是随机过程的统计平均概念，平稳随机过程的功率谱是一个确定函数；
频谱是随机过程样本的Fourier变换，对于一个随机过程而言，频谱也是一个“随机过程”。包含了幅值和相位信息。
能量谱是频谱模的平方。能量谱密度的定义适用于能量集中在一个时间窗口附近的瞬变（脉冲状信号）。平稳信号持续无限长时间。这就是在这些情况下不能使用上面定义的能量谱密度的原因。
通常所指的能量谱和能量谱密度是一个概念；功率谱和功率谱密度是一个概念，而且功率是指平均功率。
功率谱可以从两方面来定义，一个是自相关函数的傅立叶变换，另一个是时域信号傅氏变换模平方然后除以时间长度。第一种定义就是常说的维纳辛钦定理，而第二种其实从能量谱密度来的。根据parseval定理，信号傅氏变换模平方被定义为能量谱，能量谱密度在时间上平均就得到了功率谱。
平稳信号功率谱密度恰好是其自相关函数的傅氏变换。对于非平稳信号，其自相关函数的时间平均（对时间积分，随时变性消失而再次退变成一维函数）与功率谱密度仍是傅氏变换对；
当信号在\(t \in (-\infty,\infty)\)的平均功率有限时，能量是无限的；当信号在\(t \in (-\infty,\infty)\)的能量有限时，其平均功率为0。能量有限的信号称为能量信号；平均功率有限的信号称为功率信号。
1.3.1 频谱
\[ \lim_{T \to \infty}\hat{f_T}(w) \]
1.3.2 能量谱密度
\[ \lim_{T \to \infty}{|\hat{f_T}(w)|}^2 \]
1.3.3 功率谱密度
\[ \lim_{T \to \infty}\frac{{|\hat{f_T}(w)|}^2}{2T} \]
1.4 单位阶跃信号
单位阶跃函数，又称赫维赛德阶跃函数，定义如下： \[ H(x)={\begin{cases}0,&x<0\\{\frac {1}{2}},&x=0\\1,&x>0\end{cases}} \]
或 \[ H(x)={\frac {1}{2}}\left(1+\operatorname{sgn}(x)\right) \]
它是个不连续函数，其“微分”是狄拉克 δ 函数。它是一个几乎必然是零的随机变数的累积分布函数。
事实上，\(x = 0\) 的值在函数应用上并不重要，可以任意取。这个函数由奥利弗·赫维赛德提出。
H=function(x){
if(x<0) {return (0)}
if(x==0) {return (1/2)}
if(x>0) {return (1)}
}
x <- seq(-10, 10, by=0.01)
y <- sapply(x,H)
dh <- data.frame(x, y)
g <- ggplot(dh)
g1 <- g + geom_line(col='red',mapping=aes(x,y))+
xlim(-10,10)+ylim(-1,3)+
ylab("H")+ggtitle("H(x)")
print(g1)
2 单位阶跃信号的谱分析
2.1 单位直流信号频谱
\(x(t)=1\)
由于该信号 \(t \in (-\infty,\infty)\) ，所以必须乘以双边实指数衰减函数 \(e^{-a|t|}\) 以保证信号是绝对可积。
\[ \begin{aligned} F(1) &= \lim_{a \to 0} F\left(1 \cdot e^{-a|t|} \right) \\ &= \lim_{a \to 0}\int_{-\infty}^{\infty} e^{-a|t|}e^{-jwt} dt \\ &= \lim_{a \to 0}\left(\frac{2a}{a^2+w^2}\right) \\ \end{aligned} \]
故有 \[ F(1)= \begin{cases} 0, &w\neq 0\\ \infty, &w = 0 \end{cases} \]
因为 \(F(1)\) 在 \(w=0\) 处存在 \(\delta(w)\) , 其强度为 \[ \lim_{a \to 0}\left(\frac{2a}{a^2+w^2}\right)=2tg^{-1}(\frac{w}{a})|_{-\infty}^{\infty}=2\pi \]
故可以写成 \[ F(1)=2\pi\delta(w) \]
这就是广义傅里叶变换。
2.2 符号函数信号(signum)频谱
\[ x(t)=sgn(t)=\begin{cases} -1, &t<0\\ 0, &t=0\\ 1, &t>0 \end{cases} \]
通过广义傅里叶变换，有
\[ \begin{aligned} F[sgn(t)] &=\lim_{a \to 0}\int_{-\infty}^{\infty} e^{-a|t|}sgn(t)e^{-jwt} dt \\ &=\lim_{a \to 0}\int_{-\infty}^{0}-e^{(a-jw)t}dt+\lim_{a \to 0}\int_{0}^{\infty}e^{-(a+jw)t}dt \\ &=\lim_{a \to 0}\left[\frac{-1}{a-jw}+\frac{1}{a+jw}\right] \\ &=\begin{cases} -j\frac{2}{w}, &w \neq 0\\ 0, &w = 0 \end{cases} \end{aligned} \]
2.3 单位阶跃信号频谱
\[ u(t)={\frac {1}{2}}\left(1+sgn(t)\right) \]
故 \[ F[u(t)]=\begin{cases} \frac{1}{2}[2\pi\delta(w)], &w = 0 \\ \frac{1}{2}[2\pi\delta(w)-j\frac{2}{w}], &w \neq 0 \end{cases} \]
即 \[ F[u(t)]=\begin{cases} \pi\delta(w), &w = 0 \\ -j\frac{1}{w}, &w \neq 0 \end{cases} \]
2.4 单位阶跃信号能量谱密度
\[ \begin{aligned} E(w)=|F[u(t)]|^2 &=\begin{cases} |\pi\delta(w)|^2, &w = 0 \\ |-j\frac{1}{w}|^2, &w \neq 0 \end{cases} \\ &=\begin{cases} \pi^2\delta^2(w), &w = 0 \\ \frac{1}{w^2}, &w \neq 0 \end{cases} \end{aligned} \]
E=function(w){
if(w==0) {return (Inf)}
else{return (1/w^2)}
}
w <- seq(0, 1, by=0.001)
y <- sapply(w,E)
dh <- data.frame(w, y)
g <- ggplot(dh)
g1 <- g + geom_line(col='red',mapping=aes(w,y))+
xlim(0,1)+ylim(0,100)+
ylab("E")+ggtitle("E(w)")
print(g1)
2.5 单位阶跃信号频率谱密度
\[ \mbox{频率谱密度}=\lim_{T \to \infty}\frac{\mbox{能量谱密度}}{T} \]
3 谱估计
功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.
- 经典谱估计方法
- 直接法（周期图法）：先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计
- 间接法：先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计
- 分段平均周期图法
- 加窗平均周期图法
- 现代谱估计方法
- 参数模型谱估计
- AR模型
- MA模型
- ARMA模型等
- 非参数模型谱估计
- 最小方差法
- MUSIC法等
- 参数模型谱估计
直接法（周期图法），Welth法，Barlett法， 根据定义通过自相关函数求的方法，Cooley-Turkey法
3.1 各种谱估计方法的实现
3.1.1 直接法（周期图法）
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
\(\{X_n\}\)的傅里叶变化为： \[ x(w_n) = \sum_{k=0}^{N-1} X_k e^{-i k w_n} \qquad n = 0,\dots,N-1 \] 其中 \(w_n=\frac{2\pi n}{N}\)
\[ \hat{f}(w_n)=\frac{x(w_n)\cdot \bar{x}(w_n)}{N} \]
由于 \(\hat{f}(w_n)\) 在 \(n = 0,\dots,N-1\) 是对称的，故可去掉后一半。又由于周期图在 0 频率反应的是样本均值，不是周期，在分析问题时通常不予考虑。
direct=function(x,Fs=1){
N=length(x)
fft_x=fft(x)
fft_x_conj=Conj(fft_x)
y=as.numeric(fft_x_conj*fft_x/N)[1:(N/2+1)][-1]
w=((1:N-1)*Fs/N)[1:(N/2+1)][-1]
dh <- data.frame(w, y)
g <- ggplot(dh)
g1 <- g + geom_line(col='red',mapping=aes(w,y))+
scale_y_log10()+
ylab("f")+ggtitle("直接周期图法")
return(g1)
}3.1.2 间接法
间接法先由序列 \(X(n)\) 估计出自相关函数 \(\gamma_k\) ，然后对 \(\gamma_k\) 进行傅立叶变换，便得到 \(X(n)\) 的功率谱估计。
\[ \hat{f}(w)=-\hat{\gamma_0} + 2\sum_{k=0}^{N-1}\hat{\gamma_k}\cos{wk} \]
该估计不是一致的
gamma=function(z){
Ez=mean(z)
n=length(z)
Gamma=array(NA,n-1)
sum=0
for(t in 1:n)
sum=sum+(z[t]-Ez)*(z[t]-Ez)
Gamma0=1/n*sum
for(k in 1:(n-1)){
sum=0
for(t in 1:(n-k))
sum=sum+(z[t]-Ez)*(z[t+k]-Ez)
Gamma[k]=1/n*sum
}
return(c(Gamma0,Gamma))
}
indirect=function(a){
N=length(a)
Fs=1
Gamma=gamma(a)
fft_Gamma=fft(Gamma)
f=(-Gamma[1]+2*Re(fft_Gamma))[1:(N/2+1)][-1]
w=((1:N-1)*Fs/N)[1:(N/2+1)][-1]
y=f
dh <- data.frame(w, y)
g <- ggplot(dh)
g1 <- g + geom_line(col='red',mapping=aes(w,y))+
scale_y_log10()+
ylab("f")+ggtitle("间接法")
return(g1)
}3.2各种模型的谱估计
3.2.1 高斯白噪声谱估计
生成高斯白噪声
set.seed("1234")
n=200
a=rnorm(n)
qplot(1:n,a,geom = "point",col = "red")
高斯白噪声谱估计
#直接周期图法
grid.arrange(direct(a), indirect(a), nrow=1, ncol=2)
#基础包里的几种算法
par(mfcol=c(1,3))
spectrum(a)
spectrum(a, spans = 3)
spectrum(a,method = "ar")
par(mfcol=c(1,1))3.2.2 阶跃信号谱估计
取样本为50个0,50个1，频率为1，估计阶跃样本的谱
生成阶跃信号
#测试样本
N=200
Fs=1
x = c(rep(0,N/2),rep(1,N/2))
a=data.frame(1:N,x)
g0 <- ggplot(a,mapping = aes(a[,1],a[,2]))+geom_point(col='red')+
ylab("x")+xlab("n")+ggtitle("测试样本")
print(g0)
谱估计
grid.arrange(direct(x), indirect(x), nrow=1, ncol=2)
因为有负数产生导致没画出来，为啥间接法会有负数产生,并且都是很小的负数（\(10^{-16}\) 量级），是机器误差吗.跟自带函数差别这么大。。
par(mfcol=c(1,3))
spectrum(x)
spectrum(x, spans = 3)
spectrum(x,method = "ar")
par(mfcol=c(1,1))3.2.3 加有白噪声的阶跃信号谱估计
取样本为50个0,50个10，然后加入白噪声，频率为1，进行谱估计。看是不是上一部分是不是机器误差引起的
生成信号
#白噪声
set.seed("1234")
N=200
a=rnorm(N)
#阶跃信号
Fs=1
x = c(rep(0,N/2),rep(10,N/2))
x_a=x+a
a=data.frame(1:N,x_a)
g0 <- ggplot(a,mapping = aes(a[,1],a[,2]))+geom_point(col='red')+
ylab("x")+xlab("n")+ggtitle("测试样本")
print(g0)
谱估计
grid.arrange(direct(x_a), indirect(x_a), nrow=1, ncol=2)
还是一模一样，应该是计算机机器误差引起的吧
par(mfcol=c(1,3))
spectrum(x_a)
spectrum(x_a, spans = 3)
spectrum(x_a,method = "ar")
par(mfcol=c(1,1))3.2.4 加有白噪声的尖角信号谱估计
取样本为50个0,50个10，然后加入白噪声，频率为1，进行谱估计。看是不是上一部分是不是机器误差引起的
生成信号
#白噪声
set.seed("1234")
N=200
a=rnorm(N)
#尖角信号
Fs=1
x = c(0.1*(1:100),10-0.1*(1:100))
x_a=x+a
a=data.frame(1:N,x_a)
g0 <- ggplot(a,mapping = aes(a[,1],a[,2]))+geom_point(col='red')+
ylab("x")+xlab("n")+ggtitle("测试样本")
print(g0)
谱估计
grid.arrange(direct(x_a), indirect(x_a), nrow=1, ncol=2)
par(mfcol=c(1,3))
spectrum(x_a)
spectrum(x_a, spans = 3)
spectrum(x_a,method = "ar")
par(mfcol=c(1,1))