阶跃函数的谱（2）

1 library

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

2 谱估计方法

2.1 直接法

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) 
}

3 阶跃信号谱估计

3.1 阶跃大小对谱的影响

3.1.1 样本

N=200
Fs=1
x1 = c(rep(0,N/2),rep(1,N/2))
x2 = c(rep(0,N/2),rep(2,N/2))
a1=data.frame(No=1:N,X=x1,sample='1')
a2=data.frame(No=1:N,X=x2,sample='2')

g1 <- ggplot(a1)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ylim(0,2.1)+ggtitle("测试样本1") 
g2 <- ggplot(a2)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ylim(0,2.1)+ggtitle("测试样本2")
grid.arrange(g1, g2, nrow=1, ncol=2)

3.1.2 谱估计

b1=cbind(direct(x1)$xy,sample="1")
b2=cbind(direct(x2)$xy,sample="2")
b=rbind(b1,b2)
ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")

放大一点

ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+ylab("f")+xlab("w")+ggtitle("谱估计")+
  xlim(0,0.1)

可以看出，阶跃大小变成2倍，谱估计变成4倍。但这是因为，样本比较特殊，是从0-1,0-2的跃迁，都是从0开始。

\(\{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} \]

由公式知，当 \({X_k'}=2{X_k}\) 时，\(x(w_n)'=2x(w_n)\), \(\hat{f}(w_n)'=4\hat{f}(w_n)\)。即当 \(X_k\) 整体扩大或缩小相同倍数，其谱密度变换该倍数的平方

3.2 阶跃时间对谱的影响

3.2.1 样本

N=200
Fs=1
x1 = c(rep(0,N/4),rep(1,N*3/4))
x2 = c(rep(0,N/2),rep(2,N/2))
x3 = c(rep(0,N*3/4),rep(1,N/4))
a1=data.frame(No=1:N,X=x1,sample='1')
a2=data.frame(No=1:N,X=x2,sample='2')
a3=data.frame(No=1:N,X=x3,sample='3')

g1 <- ggplot(a1)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本1") 
g2 <- ggplot(a2)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本2")
g3 <- ggplot(a3)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本3")
grid.arrange(g1, g2, g3, nrow=1, ncol=3)

3.2.2 谱估计

b1=cbind(direct(x1)$xy,sample="1")
b2=cbind(direct(x2)$xy,sample="2")
b3=cbind(direct(x3)$xy,sample="3")
b=rbind(b1,b2,b3)
ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")

放大一点

ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")+xlim(0,0.1)

只出现两条曲线的原因是测试样本1与测试样本3的谱密度估计是一样的，可以看出，在时间上对称的样本的谱也是一样的

grid.arrange(direct(x1)$g1, direct(x1)$g1_logy, direct(x2)$g1, direct(x2)$g1_logy, nrow=2, ncol=2)

#样本1
direct_x1=direct(x1)
direct_x1$g1

direct_x1$g1_logy

direct_x1$xy

##         w            y
## 1   0.005 1.013295e+01
## 2   0.010 5.067726e+00
## 3   0.015 1.126625e+00
## 4   0.020 1.200085e-31
## 5   0.025 4.061191e-01
## 6   0.030 5.645651e-01
## 7   0.035 2.076133e-01
## 8   0.040 6.548162e-33
## 9   0.045 1.259246e-01
## 10  0.050 2.043173e-01
## 11  0.055 8.457485e-02
## 12  0.060 1.211279e-32
## 13  0.065 6.079369e-02
## 14  0.070 1.050719e-01
## 15  0.075 4.587431e-02
## 16  0.080 5.940585e-33
## 17  0.085 3.590459e-02
## 18  0.090 6.423760e-02
## 19  0.095 2.891520e-02
## 20  0.100 0.000000e+00
## 21  0.105 2.382711e-02
## 22  0.110 4.357550e-02
## 23  0.115 2.000889e-02
## 24  0.120 3.035631e-32
## 25  0.125 1.707107e-02
## 26  0.130 3.170045e-02
## 27  0.135 1.476284e-02
## 28  0.140 1.783921e-34
## 29  0.145 1.291679e-02
## 30  0.150 2.425920e-02
## 31  0.155 1.141770e-02
## 32  0.160 9.418950e-35
## 33  0.165 1.018416e-02
## 34  0.170 1.929584e-02
## 35  0.175 9.157350e-03
## 36  0.180 3.466611e-33
## 37  0.185 8.293908e-03
## 38  0.190 1.582590e-02
## 39  0.195 7.561299e-03
## 40  0.200 6.310887e-32
## 41  0.205 6.934727e-03
## 42  0.210 1.331008e-02
## 43  0.215 6.395033e-03
## 44  0.220 5.067984e-34
## 45  0.225 5.927222e-03
## 46  0.230 1.143293e-02
## 47  0.235 5.519424e-03
## 48  0.240 1.377520e-33
## 49  0.245 5.162147e-03
## 50  0.250 1.000000e-02
## 51  0.255 4.847729e-03
## 52  0.260 4.699269e-34
## 53  0.265 4.569931e-03
## 54  0.270 8.886257e-03
## 55  0.275 4.323635e-03
## 56  0.280 3.358822e-33
## 57  0.285 4.104608e-03
## 58  0.290 8.008394e-03
## 59  0.295 3.909331e-03
## 60  0.300 0.000000e+00
## 61  0.305 3.734861e-03
## 62  0.310 7.309276e-03
## 63  0.315 3.578719e-03
## 64  0.320 4.348255e-33
## 65  0.325 3.438812e-03
## 66  0.330 6.748760e-03
## 67  0.335 3.313362e-03
## 68  0.340 3.206924e-33
## 69  0.345 3.200853e-03
## 70  0.350 6.298081e-03
## 71  0.355 3.099993e-03
## 72  0.360 1.728013e-33
## 73  0.365 3.009670e-03
## 74  0.370 5.936314e-03
## 75  0.375 2.928932e-03
## 76  0.380 3.138528e-34
## 77  0.385 2.856962e-03
## 78  0.390 5.648080e-03
## 79  0.395 2.793054e-03
## 80  0.400 6.310887e-32
## 81  0.405 2.736606e-03
## 82  0.410 5.422029e-03
## 83  0.415 2.687100e-03
## 84  0.420 6.165151e-34
## 85  0.425 2.644095e-03
## 86  0.430 5.249820e-03
## 87  0.435 2.607216e-03
## 88  0.440 2.189096e-33
## 89  0.445 2.576150e-03
## 90  0.450 5.125428e-03
## 91  0.455 2.550638e-03
## 92  0.460 1.323684e-34
## 93  0.465 2.530471e-03
## 94  0.470 5.044678e-03
## 95  0.475 2.515485e-03
## 96  0.480 8.655045e-35
## 97  0.485 2.505560e-03
## 98  0.490 5.004938e-03
## 99  0.495 2.500617e-03
## 100 0.500 0.000000e+00

#样本2
direct_x2=direct(x2)
direct_x2$g1

direct_x2$g1_logy

direct_x2$xy

##         w            y
## 1   0.005 8.106361e+01
## 2   0.010 2.251335e-31
## 3   0.015 9.012997e+00
## 4   0.020 7.543482e-31
## 5   0.025 3.248953e+00
## 6   0.030 5.581602e-32
## 7   0.035 1.660906e+00
## 8   0.040 1.074881e-31
## 9   0.045 1.007396e+00
## 10  0.050 0.000000e+00
## 11  0.055 6.765988e-01
## 12  0.060 7.553735e-32
## 13  0.065 4.863495e-01
## 14  0.070 7.838380e-31
## 15  0.075 3.669944e-01
## 16  0.080 7.366832e-33
## 17  0.085 2.872367e-01
## 18  0.090 1.962086e-32
## 19  0.095 2.313216e-01
## 20  0.100 0.000000e+00
## 21  0.105 1.906169e-01
## 22  0.110 2.588798e-33
## 23  0.115 1.600711e-01
## 24  0.120 1.328459e-31
## 25  0.125 1.365685e-01
## 26  0.130 2.312292e-32
## 27  0.135 1.181027e-01
## 28  0.140 8.450300e-33
## 29  0.145 1.033344e-01
## 30  0.150 0.000000e+00
## 31  0.155 9.134159e-02
## 32  0.160 4.302103e-33
## 33  0.165 8.147327e-02
## 34  0.170 5.235139e-32
## 35  0.175 7.325880e-02
## 36  0.180 8.575314e-33
## 37  0.185 6.635126e-02
## 38  0.190 2.564426e-32
## 39  0.195 6.049039e-02
## 40  0.200 1.009742e-30
## 41  0.205 5.547782e-02
## 42  0.210 2.943210e-34
## 43  0.215 5.116027e-02
## 44  0.220 2.490322e-32
## 45  0.225 4.741777e-02
## 46  0.230 1.720999e-32
## 47  0.235 4.415539e-02
## 48  0.240 4.061026e-33
## 49  0.245 4.129718e-02
## 50  0.250 0.000000e+00
## 51  0.255 3.878183e-02
## 52  0.260 1.232595e-33
## 53  0.265 3.655945e-02
## 54  0.270 2.378909e-32
## 55  0.275 3.458908e-02
## 56  0.280 6.118640e-33
## 57  0.285 3.283686e-02
## 58  0.290 1.990658e-32
## 59  0.295 3.127465e-02
## 60  0.300 0.000000e+00
## 61  0.305 2.987889e-02
## 62  0.310 5.348726e-33
## 63  0.315 2.862975e-02
## 64  0.320 9.089459e-32
## 65  0.325 2.751050e-02
## 66  0.330 2.239755e-32
## 67  0.335 2.650689e-02
## 68  0.340 2.810972e-33
## 69  0.345 2.560683e-02
## 70  0.350 0.000000e+00
## 71  0.355 2.479994e-02
## 72  0.360 3.366179e-33
## 73  0.365 2.407736e-02
## 74  0.370 3.072039e-33
## 75  0.375 2.343146e-02
## 76  0.380 1.976092e-33
## 77  0.385 2.285569e-02
## 78  0.390 3.178546e-33
## 79  0.395 2.234443e-02
## 80  0.400 1.009742e-30
## 81  0.405 2.189285e-02
## 82  0.410 1.957825e-32
## 83  0.415 2.149680e-02
## 84  0.420 1.553474e-32
## 85  0.425 2.115276e-02
## 86  0.430 6.434707e-33
## 87  0.435 2.085773e-02
## 88  0.440 6.278428e-34
## 89  0.445 2.060920e-02
## 90  0.450 0.000000e+00
## 91  0.455 2.040511e-02
## 92  0.460 6.354572e-34
## 93  0.465 2.024377e-02
## 94  0.470 1.532733e-32
## 95  0.475 2.012388e-02
## 96  0.480 4.036195e-32
## 97  0.485 2.004448e-02
## 98  0.490 2.809911e-35
## 99  0.495 2.000494e-02
## 100 0.500 0.000000e+00

3.3 阶跃时间对谱的影响2

已经知道阶跃时间对称，谱相同，下面调整阶跃的时间

3.3.1 样本

N=200
Fs=1

x1 = c(rep(0,N/5),rep(1,N-N/5))
x2 = c(rep(0,N/4),rep(1,N-N/4))
x3 = c(rep(0,round(N/3)),rep(1,N-round(N/3)))
x4 = c(rep(0,N/2),rep(1,N-N/2))
a1=data.frame(No=1:N,X=x1,sample='1')
a2=data.frame(No=1:N,X=x2,sample='2')
a3=data.frame(No=1:N,X=x3,sample='3')
a4=data.frame(No=1:N,X=x4,sample='4')

g1 <- ggplot(a1)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本1") 
g2 <- ggplot(a2)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本2")
g3 <- ggplot(a3)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本3")
g4 <- ggplot(a4)+geom_point(mapping = aes(No,X,col=sample))+
  ylab("x")+xlab("n")+ggtitle("测试样本4")
grid.arrange(g1, g2, g3, g4, nrow=2, ncol=2)

3.3.2 谱估计

b1=cbind(direct(x1)$xy,sample="1")
b2=cbind(direct(x2)$xy,sample="2")
b3=cbind(direct(x3)$xy,sample="3")
b4=cbind(direct(x4)$xy,sample="4")
b=rbind(b1,b2,b3,b4)
ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")

放大一点

ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")+xlim(0,0.1)+ 
  annotate('text', x = 0, y = 20, label = 't=N/2') +
  annotate('text', x = 0, y = 15, label = 't=N/3') +
  annotate('text', x = 0, y = 10, label = 't=N/4') +
  annotate('text', x = 0, y = 7, label = 't=N/5') 

4 思考

阶跃发生时间的不同，对谱密度估计会有影响，把所有发生阶跃的时间的谱密度的峰值连起来好像挺平滑的。

xx=function(n,A=1,N=200,Fs=1){
  x = c(rep(0,round(N/n)),rep(A,N-round(N/n)))
  return(x)
}

b=data.frame(stringsAsFactors = F)
for(i in 2:50){
  temp=cbind(direct(xx(i))$xy,sample=as.character(i))
  b=rbind(b,temp)
}

ggplot(b)+geom_line(mapping = aes(w,y,col=sample))+
  ylab("f")+xlab("w")+ggtitle("谱估计")+xlim(0,0.1)

挺漂亮的。。