1 library

library(ggplot2)
library(stats)
library(base)
library(gridExtra)

2 \(\delta\)函数定义

狄拉克δ函数可大概认为是实直线上的一个函数，它在原点以外的所有点函数值为0，只在原点为无穷：

\[ \delta (x)={\begin{cases}+\infty ,&x=0\\0,&x\neq 0\end{cases}} \]

并且满足约束条件

\[ \int _{-\infty }^{\infty }\delta (x)\,dx=1. \]

3 傅里叶变换

\[ F[\delta(t)]=\int_{-\infty}^{\infty}\delta(t)e^{-jwt}dt=1 \]

4 频谱

所以\(\delta(t)\)的频谱\(d(w)\)为1

y=function(x){
  return(1)
}
s=ggplot()+stat_function(mapping = aes(x=-1:1), fun=y)+ylim(0,2)+xlab("w")+ylab("d(w)")
print(s)

5 功率谱

\[ f(w)=\lim_{T \to \infty}\frac{|d(w)|^2}{T}=\lim_{T \to \infty}\frac{1}{T} \]

6 功率谱估计

6.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)
  gg <- ggplot(dh)
  gg1 <- gg + geom_line(col='red',mapping=aes(w,y))+
    #scale_y_log10()+
    ylab("f")+ggtitle("谱密度（直接周期图法）")
  gg1_logy <- gg + geom_line(col='red',mapping=aes(w,y))+
    scale_y_log10()+
    ylab("f")+ggtitle("谱密度（直接周期图法）")
  res=list(g=gg,g1=gg1,g1_logy=gg1_logy,xy=dh)
  return(res) 
}

6.2 样本

#测试样本
N=20
Fs=1

xx=function(n=1/2,A=1,N=200){
  x = rep(0,N)
  x[round(N*n)]=A
  return(x)
}

a=data.frame(1:N,xx(1/2,1,20))

g0 <- ggplot(a,mapping = aes(a[,1],a[,2]))+geom_point(col='red')+
  ylab("x")+xlab("n")+ggtitle("测试样本")
print(g0)

a=direct(xx(1/2,1,200))
print(a$g1)

6.2.1 突变大小，样本数量一定，改变突变位置

b1=data.frame(stringsAsFactors = F)
for(i in 1:20){
  temp=cbind(direct(xx(i/20,1,20))$xy,sample=as.character(i))
  b1=rbind(b1,temp)
}
ggplot(b1)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计（改变突变位置）")

6.2.2 突变位置，样本数量一定，改变突变大小（A=1，2，3…10）

b2=data.frame(stringsAsFactors = F)
for(i in 1:10){
  temp=cbind(direct(xx(1/2,i,20))$xy,sample=as.character(i))
  b2=rbind(b2,temp)
}
ggplot(b2)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计（改变突变大小）")

6.2.3 突变位置，大小一定，改变样本数量(N=10,20,30…100)

b3=data.frame(stringsAsFactors = F)
for(i in (1:10)*10){
  temp=cbind(direct(xx(1/2,1,i))$xy,sample=as.character(i))
  b3=rbind(b3,temp)
}
ggplot(b3)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计（改变样本数量）")

6.3 分析

对于有限的冲击函数

  • 突变位置的改变不会影响谱估计大小
  • 突变大小的改变会影响谱估计大小
  • 样本数量的改变会影响谱估计大小。但 \(f(w)*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}\)

因为 \[ X_k={\begin{cases}0 ,&k=0,1,\dots,N-1 \& k\neq j \\A ,&k=j \end{cases}} \]

所以 \[ \begin{aligned} x(w_n) &= \sum_{k=0}^{N-1} X_k e^{-i k w_n} \qquad n = 0,\dots,N-1 \\ &= X_j e^{-i j w_n} \end{aligned} \]

故 \[ \begin{aligned} \hat{f}(w_n) &=\frac{x(w_n)\cdot \bar{x}(w_n)}{N} \\ &=\frac{X_j^2}{N} \\ &=\frac{A^2}{N} \end{aligned} \]

这个式子可以解释刚才由谱估计分析得到的三个结论