精华内容
下载资源
问答
  • 应用循环自相关函数和快速谱峭度相结合的方法,对滚动轴承早期故障诊断进行分析研究。首先利用谱峭度方法确定滚动轴承振动信号的最佳带通滤波器,然后利用循环自相关函数对滤波后的信号进行解调,提取出滚动轴承故障...
  • 循环自相关函数,matlab程序

    热门讨论 2010-05-07 09:33:48
    自己写的关于循环自相关函数的快速计算方法,不同的调制方式只要更改x_t表达式即可
  • OFDM信号的循环自相关函数原理详解及matlab实现

    千次阅读 热门讨论 2018-03-05 18:32:12
    对接收信号的自相关函数进行FFT变换后得到OFDM信号的循环自相关函数 因此通过上述分析我们可以通过求解循环自相关函数求解出OFDM 信号的有效数据长度 和符号总长度 ,因此可以间接的求出OFDM信号的子载波间隔 和...

    对接收信号的自相关函数进行FFT变换后得到OFDM信号的循环自相关函数

        因此通过上述分析我们可以通过求解循环自相关函数求解出OFDM 信号的有效数据长度 和符号总长度 ,因此可以间接的求出OFDM信号的子载波间隔 和循环前缀的长度 。

        由循环自相关函数的特点可知,通过估计循环频率a=0截面上次峰值的位置可以估计出信号的有效数据长度 ,通过搜索 截面上的峰值间隔可以估计出信号的符号总长度。但是在求解过程中不需要直接去求解循环自相关,增加不必要的运算量。只需要求解 和a=0两个截面即可。

        首先通过可变延时自相关的方法求解出a=0的截面估计出有效数据长度 然后根据令可变延时 ,对 做FFT变换估计出符号总长度 ,进而可以算出子载波间隔和循环前缀的长度。

    给出matlab代码如下:

     function [Tu,Ts] = auto_xcorr(data, P, xcorr_len, N,t,K)
     %*************************************************************************
     %功能:计算发送数据或接收数据自相关
     %data:输入数据
     %P:循环周期
     %N:OFDM符号数
     %xcorr_len:自相关长度,以OFDM符号为单位
     %K=1时绘制图形  
     %Rx自相关函数
     %*************************************************************************
     kk = xcorr_len+2 ;
     xcorr_len = xcorr_len*P;
     Rx = zeros(xcorr_len, 2*P);
     data1=[zeros(1,P),data,zeros(1,P)];
     for i = 1 : N-kk+1
        for tao = 1-P : P
            for n = 1 : xcorr_len
                Rx(n,tao+P) = Rx(n,tao+P) + data1(P*(i-1)+n+P)*conj(data1(P*(i-1)+n+tao+P));
            end
        end
     end
    
     Rx = Rx./(N-kk+1);
     R_tao =zeros(1,length(Rx(1,:)));
     for i=1:length(Rx(:,1))
        R_tao= R_tao+Rx(i,:);
     end
    
     R_tao= abs(R_tao)/length(Rx(:,1));
     
     if K==1
         figure
         stem ((1-P)*t :t:( P)*t ,R_tao);   %P*t是将由采样点数表示的有效数据长度变为由时间表示的
         xlabel('delay/s');
         ylabel('amplitute');
         title('短波可变延时自相关')
         figure
         stem ((1-P) :( P) ,R_tao);   %P*t是将由采样点数表示的有效数据长度变为由时间表示的
         xlabel('delay/s');
         ylabel('amplitute');
         title('基带可变延时自相关')
     end
     %********************************************
     %取有效数据长度
     %********************************************
     abc=sort( R_tao(1,1:end/2));    %按升序排序
     Z=length(abc);
     [as1 ,Tu1] =find( R_tao(1,1:end/2)==abc(1,Z-1));%取峰值高度为第二高的峰值对应坐标
     Tu=(length(R_tao)/2-Tu1)*t;       %为信号的有效数据长度
    
     %********************************************
     %取数据的总长度
     %********************************************
     
     sigma0 = 6;   %观测数据的长度因子
     L_fft=sigma0*P;     %观测数据的长度,即fft的长度
     RR = zeros(xcorr_len, L_fft);
     data111=[zeros(1,P),data,zeros(1,P)];   %将观测的数据长度延长,避免求自相关时超出长度
     for i = 1 : N-4 
        for k=1:L_fft
                RR(i,k) =data111(P*(i-1)+k)*conj(data111(P*(i-1)+k+length(R_tao)/2-Tu1));
        end
     end
    
     RRR_tao =zeros(1,length(RR(1,:)));
     for i=1:length(RR(:,1))
        RRR_tao= RRR_tao+RR(i,:);
     end
    
     RRR_tao= abs(RRR_tao)/length(RR(:,1));
     RRR_T=ifft(RRR_tao);
     RRR_fft= abs(RRR_T(1,1:20*sigma0));
     if K==1
         figure
         LL=t:t:20*sigma0*t;
         stem(LL,RRR_fft);
         title('短波固定延时自相关的IFFT')
         figure
         stem(RRR_fft);
         title('基带固定延时自相关的IFFT')
     end
    
    a=sort(RRR_fft);    %按升序排序
    Z=length(a);
    [a1 ,T1] =find(RRR_fft==a(1,Z));%取峰值高度为第一高的峰值对应坐标
    [a2 ,T2] =find(RRR_fft==a(1,Z-1));%取峰值高度为第二高的峰值对应坐标
    [a3 ,T3] =find(RRR_fft==a(1,Z-2));%取峰值高度为第三高的峰值对应坐标
    Ts1 = abs(T2-T1);
    Ts2 = abs(T3-T2);
    Ts3 = abs(T3-T1);
    Tss = [Ts1 Ts2 Ts3];
    TT = sort(Tss);  %按升序排序取间隔最小值
    Ts= L_fft/TT(1)*t;       %为信号的有效数据长度
    给出参考的结果如下:


    展开全文
  • 说明:接上一节循环自相关函数和谱相关密度(二)——实信号、复信号模型下的BPSK信号循环谱推导 8 QPSK信号谱相关密度函数 8.1 实信号模型 QPSK实信号表达式可以写为 r(t)=yI(t)−yQ(t)+n(t)r(t) = {y_I}(t) - {y_Q...

    说明:接上一节循环自相关函数和谱相关密度(三)——实信号、复信号模型下的BPSK信号循环谱MATLAB仿真结果及代码

    8 QPSK信号谱相关密度函数

    8.1 实信号模型

    QPSK实信号表达式可以写为
    r(t)=yI(t)yQ(t)+n(t)r(t) = {y_I}(t) - {y_Q}(t) + n(t)

    =c(t)pI(t)s(t)pQ(t)+n(t)= c(t){p_I}(t) - s(t){p_Q}(t) + n(t)

    =n=c(nT)q(tnTt0)cos(2πf0t+θ)= \sum\limits_{n = - \infty }^\infty {c(nT)q(t - nT - {t_0})} \cos (2\pi {f_0}t + \theta )

    +n=s(nT)q(tnTt0)sin(2πf0t+θ) + n(t)+ \sum\limits_{n = - \infty }^\infty {s(nT)q(t - nT - {t_0})} \sin (2\pi {f_0}t + \theta ){\text{ + }}n(t) (54)

    其中,t0{t_0}为起始时间,TT为符号速率,c(n)c(n)s(n)s(n)为基带符号序列,f0{f_0}为载波频率,θ\theta为初始相位,n(t)n(t)为高斯白噪声,q(t)q(t)为矩形脉冲,且

    c(t)=n=a(nT)q(tnTt0)c(t) = \sum\limits_{n = - \infty }^\infty {a(nT)q(t - nT - {t_0})}

    =q(tt0)na(t)δ(tnT)= q(t - {t_0}) \otimes \sum\limits_n {a(t)\delta (t - nT)}

    =q(tt0)a^(t)= q(t - {t_0}) \otimes \hat a(t) (55)

    s(t)=n=b(nT)q(tnTt0)s(t) = \sum\limits_{n = - \infty }^\infty {b(nT)q(t - nT - {t_0})}

    =q(tt0)nb(t)δ(tnT)= q(t - {t_0}) \otimes \sum\limits_n {b(t)\delta (t - nT)}

    =q(tt0)b^(t)= q(t - {t_0}) \otimes \hat b(t) (56)

    pI(t)=cos(2πf0t+θ){p_I}(t) = \cos (2\pi {f_0}t + \theta ) (57)

    pQ(t)=sin(2πf0t+θ){p_Q}(t) = \sin (2\pi {f_0}t + \theta ) (58)

    考虑pQ(t){p_Q}(t)的二次变换

    vQτ(t)=pQ(t + τ/2)pQ(tτ/2){v_{Q\tau }}(t) = {p_Q}(t{\text{ + }}\tau /2)p_Q^*(t - \tau /2)

    =14(ej2πf0τ+ej2πf0τej(4πf0t+2θ)ej(4πf0t+2θ))= \frac{1}{4}({e^{j2\pi {f_0}\tau }} + {e^{ - j2\pi {f_0}\tau }} - {e^{j(4\pi {f_0}t + 2\theta )}} - {e^{ - j(4\pi {f_0}t + 2\theta )}}) (59)

    其Fourier级数系数为

    vQτ(t)ej2παtt=14ej2πf0τej2παtt+14ej2πf0τej2παtt{\left\langle {{v_{Q\tau }}(t){e^{ - j2\pi \alpha t}}} \right\rangle _t} = \frac{1}{4}{e^{j2\pi {f_0}\tau }}{\left\langle {{e^{ - j2\pi \alpha t}}} \right\rangle _t} + \frac{1}{4}{e^{ - j2\pi {f_0}\tau }}{\left\langle {{e^{ - j2\pi \alpha t}}} \right\rangle _t}

    14ej2θej2π(α+2f0)tt14ej2θej2π(α2f0)tt- \frac{1}{4}{e^{ - j2\theta }}{\left\langle {{e^{ - j2\pi (\alpha + 2{f_0})t}}} \right\rangle _t} - \frac{1}{4}{e^{j2\theta }}{\left\langle {{e^{ - j2\pi (\alpha - 2{f_0})t}}} \right\rangle _t} (60)

    pQ(t){p_Q}(t)的循环自相关函数和谱相关密度函数为

    RpQα(τ)={14e±j2θα=±2f012cos(2πf0τ)α=00 otherwise R_{p_{Q}}^{\alpha}(\tau)=\left\{\begin{array}{cc}-\frac{1}{4} e^{\pm j 2 \theta} & \alpha=\pm 2 f_{0} \\ \frac{1}{2} \cos \left(2 \pi f_{0} \tau\right) & \alpha=0 \\ 0 & \text { otherwise }\end{array}\right.(61)

    SpQα(f)={14e±j2θδ(f)α=±2f014[δ(f+f0)+δ(ff0)]α=00 otherwise S_{p_{Q}}^{\alpha}(f)=\left\{\begin{array}{cc}-\frac{1}{4} e^{\pm j 2 \theta} \delta(f) & \alpha=\pm 2 f_{0} \\ \frac{1}{4}\left[\delta\left(f+f_{0}\right)+\delta\left(f-f_{0}\right)\right] & \alpha=0 \\ 0 & \text { otherwise }\end{array}\right.$(62)

    由(12)、(13)得yQ(t){y_Q}(t)的循环自相关函数和谱相关密度函数为

    RyQα(τ)=βRpQβ(τ)Rsαβ(τ)R_{{y_Q}}^\alpha (\tau ) = \sum\limits_\beta {R_{{p_Q}}^\beta (\tau )R_s^{\alpha - \beta }(\tau )}

    =14ej2θRsα2f0(τ)14ej2θRsα+2f0(τ)+12cos(2πf0τ)Rsα(τ)= - \frac{1}{4}{e^{j2\theta }}R_s^{\alpha - 2{f_0}}(\tau ) - \frac{1}{4}{e^{ - j2\theta }}R_s^{\alpha + 2{f_0}}(\tau ) + \frac{1}{2}\cos (2\pi {f_0}\tau )R_s^\alpha (\tau ) (63)

    SyQα(f)=βSpQβ(f)Ssαβ(f)S_{{y_Q}}^\alpha (f) = \sum\limits_\beta {S_{{p_Q}}^\beta (f) \otimes S_s^{\alpha - \beta }(f)}

    =14[Ssα(f+f0)+Ssα(ff0)ej2θSsα2f0(f)ej2θSsα+2f0(f)]= \frac{1}{4}\left[ {S_s^\alpha (f + {f_0}) + S_s^\alpha (f - {f_0}) - {e^{j2\theta }}S_s^{\alpha - 2{f_0}}(f) - {e^{ - j2\theta }}S_s^{\alpha + 2{f_0}}(f)} \right] (64)

    由(35)可得yI(t){y_I}(t)的谱相关密度函数为
    SyIα(f)=βSpIβ(f)Scαβ(f)S_{{y_I}}^\alpha (f) = \sum\limits_\beta {S_{{p_I}}^\beta (f) \otimes S_c^{\alpha - \beta }(f)}

    =14[Scα(f+f0)+Scα(ff0)+ej2θScα2f0(f)+ej2θScα+2f0(f)]= \frac{1}{4}\left[ {S_c^\alpha (f + {f_0}) + S_c^\alpha (f - {f_0}) + {e^{j2\theta }}S_c^{\alpha - 2{f_0}}(f) + {e^{ - j2\theta }}S_c^{\alpha + 2{f_0}}(f)} \right] (65)

    同(29),可得

    Scα(f)=1TQ(f+α/2)Q(fα/2)ej2παt0S~aα(f)S_c^\alpha (f) = \frac{1}{T}Q(f + \alpha /2){Q^*}(f - \alpha /2){e^{ - j2\pi \alpha {t_0}}}\tilde S_a^\alpha (f) (66)

    Ssα(f)=1TQ(f+α/2)Q(fα/2)ej2παt0S~bα(f)S_s^\alpha (f) = \frac{1}{T}Q(f + \alpha /2){Q^*}(f - \alpha /2){e^{ - j2\pi \alpha {t_0}}}\tilde S_b^\alpha (f) (67)

    将(66)、(67)分别代入式(65)、(64),得y(t)=yI(t)yQ(t)y(t) = {y_I}(t) - {y_Q}(t)的谱相关密度函数为

    Syα(f)=12T[Q(f+f0+α/2)Q(f+f0α/2)S~aα(f+f0)S_y^\alpha (f) = \frac{1}{{2T}}[Q(f + {f_0} + \alpha /2){Q^*}(f + {f_0} - \alpha /2)\tilde S_a^\alpha (f + {f_0})

    +Q(ff0+α/2)Q(ff0α/2)S~aα(ff0)]ej2παt0+ Q(f - {f_0} + \alpha /2){Q^*}(f - {f_0} - \alpha /2)\tilde S_a^\alpha (f - {f_0})]{e^{ - j2\pi \alpha {t_0}}} (68)

    其中,假设同相和正交分量是平衡的,即

    R~a(0)=R~b(0)=1{\tilde R_a}(0) = {\tilde R_b}(0) = 1 (69)

    S~aα(f)=S~bα(f)\tilde S_a^\alpha (f) = \tilde S_b^\alpha (f) (70)

    由(38),且对于高斯白噪声n(t)n(t),当且仅当α=0\alpha = 0时,其谱相关密度函数不为零,则QPSK实信号的谱相关密度函数为

    Srα(f)={Syα(f)+Snα(f),α=0Syα(f),α0S_{r}^{\alpha}(f)=\left\{\begin{array}{cc}S_{y}^{\alpha}(f)+S_{n}^{\alpha}(f), & \alpha=0 \\ S_{y}^{\alpha}(f), & \alpha \neq 0\end{array}\right.(71)

    同样,噪声n(t)n(t)只影响循环频率为零时的截面。

    8.2 复信号模型

    QPS复信号表达式可以写为
    r(t)=y(t)+n(t)r(t) = y(t) + n(t)

    =d(t)p(t)+n(t)= d(t)p(t) + n(t)

    =[c(t)+js(t)]ej(2πf0t+θ)+n(t)= [c(t) + js(t)]{e^{j(2\pi {f_0}t + \theta )}} + n(t) (72)

    其中,t0{t_0}为起始时间,TT为符号速率,c(n)c(n)s(n)s(n)为基带符号序列,f0{f_0}为载波频率,θ\theta为初始相位,n(t)n(t)为高斯白噪声,q(t)q(t)为矩形脉冲,且

    d(t)=c(t)+js(t)d(t) = c(t) + js(t) (73)

    p(t)=pI(t)+jpQ(t) = ej(2πf0t+θ)p(t) = {p_I}(t) + j{p_Q}(t){\text{ = }}{e^{j(2\pi {f_0}t + \theta )}} (74)

    式(73)、(74)中,c(t)c(t)s(t)s(t)同8.1节。

    对于01先验等概的复基带符号序列d(n)=a(n)+jb(n)d(n) = a(n) + jb(n),其循环自相关函数为

    R~dα(kT)=limN12N+1n=NN[a(nT+kT)+jb(nT+kT)][a(nT)+jb(nT)]ej2πα(n+k/2)T\tilde R_d^\alpha (kT) = \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2N + 1}}\sum\limits_{n = - N}^N {[a(nT + kT) + jb(nT + kT)]{{[a(nT) + jb(nT)]}^*}} {e^{ - j2\pi \alpha (n + k/2)T}}

    =limN12N+1n=NN{[a(nT+kT)a(nT)+b(nT+kT)b(nT)]= \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2N + 1}}\sum\limits_{n = - N}^N {\{ [a(nT + kT)a(nT) + b(nT + kT)b(nT)]}

    +j[b(nT+kT)a(nT)a(nT+kT)b(nT)]ej2πα(n+k/2)T+ j[b(nT + kT)a(nT) - a(nT + kT)b(nT)]{e^{ - j2\pi \alpha (n + k/2)T}} (75)

    同理,假设同相和正交分量是平衡的,则当且仅当k=0k = 0α=m/T\alpha = m/T时,R~dα(kT)=2R~aα(kT)=2R~bα(kT)\tilde R_d^\alpha (kT) = 2\tilde R_a^\alpha (kT) = 2\tilde R_b^\alpha (kT),则其谱相关密度函数为

    S~dα(f)={2R~a(0)=2,α=m/T0,αm/T\tilde{S}_{d}^{\alpha}(f)=\left\{\begin{aligned} 2 \tilde{R}_{a}(0)=2, & \alpha=m / T \\ 0, & \alpha \neq m / T \end{aligned}\right.(76)

    由(29)、(46)、(47),得y(t)y(t)的谱相关密度函数为

    Syα(f)=2T[Q(ff0+α/2)Q(ff0α/2)ej2παt0S~aα(ff0)]S_y^\alpha (f) = \frac{2}{T}[Q(f - {f_0} + \alpha /2){Q^*}(f - {f_0} - \alpha /2){e^{ - j2\pi \alpha {t_0}}}\tilde S_a^\alpha (f - {f_0})] (77)

    同(71),可得复QPSK信号的谱相关密度函数为

    Srα(f)={Syα(f)+Snα(f),α=0Syα(f),α0S_{r}^{\alpha}(f)=\left\{\begin{array}{cc}S_{y}^{\alpha}(f)+S_{n}^{\alpha}(f), & \alpha=0 \\ S_{y}^{\alpha}(f), & \alpha \neq 0\end{array}\right.(78)

    同样,噪声n(t)n(t)只影响循环频率为零时的截面。

    8.3 谱分析

    8.3.1 主峰个数

    对于QPSK实信号,由(68)、(70)可知,其谱相关密度函数在α=0\alpha = 0f=±f0f = \pm {f_0}处各有一个主峰,实QPSK信号共有2个主峰。

    对于复QPSK信号,由(77)、(78)可知,其谱相关密度函数仅在f=f0f = {f_0}α=0\alpha = 0处有一个谱峰。

    8.3.2 切面特征

    在式(68)中,令f=0f = 0α=±2f0+m/T\alpha = \pm 2{f_0} + m/T

    Syα(f)=12T[Q(f0+α/2)Q(f0α/2)S~aα(f0)S_y^\alpha (f) = \frac{1}{{2T}}[Q({f_0} + \alpha /2){Q^*}({f_0} - \alpha /2)\tilde S_a^\alpha ({f_0})

    +Q(f0+α/2)Q(f0α/2)S~aα(f0)]ej2παt0+ Q( - {f_0} + \alpha /2){Q^*}( - {f_0} - \alpha /2)\tilde S_a^\alpha ( - {f_0})]{e^{ - j2\pi \alpha {t_0}}} (79)

    特别地,当α=±2f0=n/T\alpha = \pm 2{f_0} = n/T时,有

    Syα(f)=12T[Q(2f0)Q(0)+Q(0)Q(2f0)]ej2πnt0/TS_y^\alpha (f) = \frac{1}{{2T}}[Q(2{f_0}){Q^*}(0) + Q(0){Q^*}( - 2{f_0})]{e^{ - j2\pi n{t_0}/T}} (80)

    一般二倍载波频率2f02{f_0}大于符号速率1/T1/T,由(80)可以看出,当f=0f = 0α=±2f0\alpha = \pm 2{f_0}附近会出现较小峰值,但在二倍载波频率2f02{f_0}为符号速率1/T1/T的整数倍时,其谱相关密度函数不一定取到极大值,即没有明显的离散谱线存在,故此时无法利用f=0f = 0切面对实QPSK信号载频进行估计。

    f=±f0f = \pm {f_0}α=m/T\alpha = m/T

    Syα(f)={12T[Q(2f0+α/2)Q(2f0α/2)+Q(α/2)2]ej2παt0f=f012T[Q(α/2)2+Q(2f0+α/2)Q(2f0α/2)]ej2παt0f=f0S_{y}^{\alpha}(f)=\left\{\begin{array}{cc}\frac{1}{2 T}\left[Q\left(2 f_{0}+\alpha / 2\right) Q^{*}\left(2 f_{0}-\alpha / 2\right)+|Q(\alpha / 2)|^{2}\right] e^{-j 2 \pi \alpha t_{0}} & f=f_{0} \\ \frac{1}{2 T}\left[|Q(\alpha / 2)|^{2}+Q\left(-2 f_{0}+\alpha / 2\right) Q^{*}\left(-2 f_{0}-\alpha / 2\right)\right] e^{-j 2 \pi \alpha t_{0}} & f=-f_{0}\end{array}\right.(81)

    即在f=±f0f = \pm {f_0}切面,其谱相关密度函数幅度在循环频率为α=m/T\alpha = m/T即符号速率整数倍处出现峰值,且频率分辨率较小时峰值较明显,由此可估计实QPSK信号的符号速率,此外还可根据符号速率处对应的谱相关密度函数的相位来估计时延t0{t_0}

    对于复QPSK信号,对比(77)和(47)可知,复QPSK信号为复BPSK信号谱相关密度函数乘以常数2,因此二者特性相同,即在α=0\alpha = 0切面,其谱相关密度函数幅度只在f=f0f = {f_0}出现峰值,由此可估计复QPSK信号的载波频率,但此时噪声n(t)n(t)的谱相关密度函数不为零,因此利用该切面进行载频估计受噪声影响较大。在f=f0f = {f_0}切面,其谱相关密度函数幅度在循环频率为α=m/T\alpha = m/T即符号速率整数倍处出现峰值,在α=0\alpha = 0处的峰值最大,由此可估计实QPSK信号的符号速率,此外还可根据符号速率处对应的谱相关密度函数的相位来估计时延t0{t_0}

    展开全文
  • 说明:接上一节循环自相关函数和谱相关密度(四)——实信号、复信号模型下的QPSK信号循环谱推导 8.4 仿真结果 8.4.1 实QPSK信号 符号速率RB = 40,采样率Fs = 960,载波频率fc = 300,符号数N = 1000,矩形成形,...

    说明:接上一节循环自相关函数和谱相关密度(四)——实信号、复信号模型下的QPSK信号循环谱推导

    8.4 仿真结果

    8.4.1 实QPSK信号

    • 符号速率RB = 40,采样率Fs = 960,载波频率fc = 300,符号数N = 1000,矩形成形,二倍载波频率为符号速率的整数倍。
      在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述
    • 符号速率RB = 60,采样率Fs = 960,载波频率fc = 400,符号数N = 1000,矩形成形,二倍载波频率不为符号速率的整数倍。
      在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述
    • 符号速率RB = 60,采样率Fs = 960,载波频率fc = 300,符号数N = 1000,滚降系数 ,根升余弦成形,二倍载波频率为符号速率的整数倍。
      在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述在这里插入图片描述

    8.4.2 复QPSK信号

    • 符号速率RB = 60,采样率Fs = 960,载波频率fc = 400,符号数N = 1000,矩形成形。
      在这里插入图片描述在这里插入图片描述
      在这里插入图片描述在这里插入图片描述
    • 符号速率RB = 60,采样率Fs = 960,载波频率fc = 400,符号数N = 1000,根升余弦成形。
      在这里插入图片描述在这里插入图片描述
      在这里插入图片描述在这里插入图片描述

    8.5 仿真代码

    DQPSK实信号

    clear;close all;clc;
    %% 生成dqpsk信号
    RB = 60; % 符号速率
    beta = 0.8; % 滚降系数
    B = (1+beta)*RB; % 带宽
    Fs = 960; % 采样率
    fc = 300; % 载频
    symlen = 1000; % 符号数
    rc = randi([0,4-1],1,symlen);
    fd = Fs/RB; % 过采倍数
    Len = symlen*fd; % 总的采样点数
    t = 0:1/Fs:(Len-1)/Fs; % 共symlen*fd个样点,1s采样960个点,时间为6.25秒
    
    % 将绝对码变为相对码,采用数组查表的方法
    drc = zeros(1,symlen); % 参考比特第1比特设置为1
    % 星座映射关系为表7-1,P260
    % 绝对码映射到相位变化Δφ 0(1)---0°; 1(2)--- 270°; 2(3)---90°; 3(4)---180°
    % 如270表示绝对码比特1代表相位变化270°,90表示比特2相位变化90°
    dphai = [0,270,90,180];
    % 当前时刻相位φ映射到星座点 0(1)---0°; 1(2)--- 270°; 2(3)---90°; 3(4)---180°
    % 如当前时刻相位为180°,则mod(180,360)/90=2,当前相位映射为phaik(2+1)=3
    % 式中加1是因为MATLAB坐标从1算起
    phaik = [0,2,3,1]; 
    for i = 2:symlen
        drc(i) = phaik(mod(dphai(rc(i)+1)+dphai(drc(i-1)+1),360)/90+1);
    end
    
    % 将四进制数据分成同相、正交两路数据的双极性码
    % I支路在低位,Q支路在高位
    drc_i = zeros(1,symlen);
    drc_q = zeros(1,symlen);
    for i = 1:symlen
        drcb = dec2bin(drc(i), 2);
        drc_i(i) = 1-2*str2double(drcb(2));
        drc_q(i) = 1-2*str2double(drcb(1));
    end
    
    %% 成形滤波
    % 矩形成形
    for i = 1:symlen
        sbas_rect_i((i-1)*fd+1:i*fd) = drc_i(i)*ones(1,fd);
        sbas_rect_q((i-1)*fd+1:i*fd) = drc_q(i)*ones(1,fd);
    end
    
    % 根升余弦成形
    h = rcosdesign(beta,5,fd,'sqrt');
    up_drc_i = upsample(drc_i,fd);
    up_drc_q = upsample(drc_q,fd);
    sbas_rcos_i = filter(h,1,up_drc_i);
    sbas_rcos_q = filter(h,1,up_drc_q);
    
    %% 调制
    % 实信号
    % dqpsk = sbas_rect_i.*cos(2*pi*fc*t)-sbas_rect_q.*sin(2*pi*fc*t); % 矩形成形
    dqpsk = sbas_rcos_i.*cos(2*pi*fc*t)-sbas_rcos_q.*sin(2*pi*fc*t); % 根升余弦成形
    figure, plot(-Fs/2:Fs/2-1, fftshift(abs(fft(dqpsk,Fs))),'k');
    title('DQPSK实信号频谱');
    N = 4096;
    M = 128;
    [Sx, alphao, fo] = autofam(dqpsk, Fs, Fs/N*M, Fs/N); % 循环谱
    figure; mesh(alphao,fo,abs(Sx)); 
    title('DQPSK实信号SCF');xlabel('Cycle frequency(alpha)');
    ylabel('frequency(f)');zlabel('Sx(f)');grid on;
    figure; contour(fo,alphao,abs(Sx).',[0.1 0.1],'k');
    axis([-1000 1000 -1000 1000]);title('SCF');
    xlabel('Frequency(f)');ylabel('Cycle frequency(alpha)');
    
    %% f = 0 截面
    f0 = (length(fo)+1 )/2;
    figure; plot(alphao,abs(Sx(f0,:)),'k');title('DQPSK实信号SCF, f=0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');ylim([0 1.1]);
    
    %% f = f0 截面
    Nf = pow2(nextpow2(N/M)); % f 轴分辨率 1/Np
    L = Nf/4;
    P = pow2(nextpow2(4 * N/Nf));
    Na = P*L; % α 轴分辨率 1/Na
    ff0 = floor((fc+Fs/2)*Nf/Fs+1); % -480+Fs/Np*(k-1) = 300
    figure; plot(alphao,abs(Sx(ff0,:)),'k');title('DQPSK实信号SCF, f=f0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');ylim([0 1.1]);
    
    %% f = -f0 截面
    ff0 = floor((-fc+Fs/2)*Nf/Fs+1); % -480+Fs/Np*(k-1) = -300
    figure; plot(alphao,abs(Sx(ff0,:)),'k');title('DQPSK实信号SCF, f=-f0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');ylim([0 1.1]);
    
    %% alpha = 0 截面
    a0 = (length(alphao)+1)/2;
    figure; plot(fo,abs(Sx(:,a0)),'k');title('DQPSK实信号SCF, α=0');
    xlabel('Frequency(f)');ylabel('Sx(f)');grid on;
    

    DQPSK复信号

    clear;close all;clc;
    %% 生成dqpsk信号
    RB = 60; % 符号速率
    beta = 0.8; % 滚降系数
    B = (1+beta)*RB; % 带宽
    Fs = 960; % 采样率
    fc = 400; % 载频
    symlen = 1000; % 符号数
    rc = randi([0,4-1],1,symlen);
    fd = Fs/RB; % 过采倍数
    Len = symlen*fd; % 总的采样点数
    t = 0:1/Fs:(Len-1)/Fs; % 共symlen*fd个样点,1s采样960个点,时间为6.25秒
    
    % 将绝对码变为相对码,采用数组查表的方法
    drc = zeros(1,symlen); % 参考比特第1比特设置为1
    % 星座映射关系为表7-1,P260
    % 绝对码映射到相位变化Δφ 0(1)---0°; 1(2)--- 270°; 2(3)---90°; 3(4)---180°
    % 如270表示绝对码比特1代表相位变化270°,90表示比特2相位变化90°
    dphai = [0,270,90,180];
    % 当前时刻相位φ映射到星座点 0(1)---0°; 1(2)--- 270°; 2(3)---90°; 3(4)---180°
    % 如当前时刻相位为180°,则mod(180,360)/90=2,当前相位映射为phaik(2+1)=3
    % 式中加1是因为MATLAB坐标从1算起
    phaik = [0,2,3,1]; 
    for i = 2:symlen
        drc(i) = phaik(mod(dphai(rc(i)+1)+dphai(drc(i-1)+1),360)/90+1);
    end
    
    % 将四进制数据分成同相、正交两路数据的双极性码
    % I支路在低位,Q支路在高位
    drc_i = zeros(1,symlen);
    drc_q = zeros(1,symlen);
    for i = 1:symlen
        drcb = dec2bin(drc(i), 2);
        drc_i(i) = 1-2*str2double(drcb(2));
        drc_q(i) = 1-2*str2double(drcb(1));
    end
    
    %% 成形滤波
    % 矩形成形
    for i = 1:symlen
        sbas_rect_i((i-1)*fd+1:i*fd) = drc_i(i)*ones(1,fd);
        sbas_rect_q((i-1)*fd+1:i*fd) = drc_q(i)*ones(1,fd);
    end
    
    % 根升余弦成形
    h = rcosdesign(beta,5,fd,'sqrt');
    up_drc_i = upsample(drc_i,fd);
    up_drc_q = upsample(drc_q,fd);
    sbas_rcos_i = filter(h,1,up_drc_i);
    sbas_rcos_q = filter(h,1,up_drc_q);
    
    %% 调制
    % 复信号
    % dqpsk_c = (sbas_rect_i+1j*sbas_rect_q).*exp(1j*2*pi*fc*t); % 矩形成形
    dqpsk_c = (sbas_rcos_i+1j*sbas_rcos_q).*exp(1j*2*pi*fc*t); % 根升余弦成形
    figure, plot(-Fs/2:Fs/2-1, fftshift(abs(fft(dqpsk_c,Fs))),'k');
    title('DQPSK复信号频谱');
    N = 4096;
    M = 64;
    [Sx, alphao, fo] = autofam(dqpsk_c, Fs, Fs/N*M, Fs/N); % 循环谱
    figure; mesh(alphao,fo,abs(Sx)); 
    title('DQPSK复信号SCF');xlabel('Cycle frequency(alpha)');
    ylabel('frequency(f)');zlabel('Sx(f)');grid on;
    figure; contour(fo,alphao,abs(Sx).',[0.1 0.1],'k');
    axis([-1000 1000 -1000 1000]);title('SCF');
    xlabel('Frequency(f)');ylabel('Cycle frequency(alpha)');
    
    %% f = 0 截面
    f0 = (length(fo)+1 )/2;
    figure; plot(alphao,abs(Sx(f0,:)),'k');title('DQPSK复信号SCF, f=0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');
    % ylim([0 1.1]);
    
    %% f = f0 截面
    Nf = pow2(nextpow2(N/M)); % f 轴分辨率 1/Np
    L = Nf/4;
    P = pow2(nextpow2(4 * N/Nf));
    Na = P*L; % α 轴分辨率 1/Na
    ff0 = floor((fc+Fs/2)*Nf/Fs+1); % -480+Fs/Np*(k-1) = 300
    figure; plot(alphao,abs(Sx(ff0,:)),'k');title('DQPSK复信号SCF, f=f0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');
    % ylim([0 1.1]);
    
    %% f = -f0 截面
    ff0 = floor((-fc+Fs/2)*Nf/Fs+1); % -480+Fs/Np*(k-1) = -300
    figure; plot(alphao,abs(Sx(ff0,:)),'k');title('DQPSK复信号SCF, f=-f0');
    xlabel('Cycle frequency(alpha)');ylabel('Sx(f)');
    % ylim([0 1.1]);
    
    %% alpha = 0 截面
    a0 = (length(alphao)+1)/2;
    figure; plot(fo,abs(Sx(:,a0)),'k');title('DQPSK复信号SCF, α=0');
    xlabel('Frequency(f)');ylabel('Sx(f)');grid on;
    

    FAM法计算循环谱子函数(未仔细研究推导该算法,计算速度挺快)

     function [Sx,alphao,fo] = autofam(x,fs,df,dalpha)
    %************************************************
    % x为输入信号向量
    % fs为采样率
    % df为频率分辨率
    % dalpha为循环频率分辨率
    % 注意:需满足df>>dalpha才能得到较高的效果,fs/dalpha为进行采样计算功率谱的码元个数
    %%fam算法比ssca(autossca)算法快,但效果一样.都是求循环平稳与谱
    %************************************************
    % if nargin > 4 | nargin < 4
    %     error('Wrong number of arguments.');
    % end   
    % fs=1280;
    % fc=160;
    % N=8192;
    % PW=N/fs;
    % signal_type=5;
    % switch signal_type
    %     case 1; input_signal = sin(2*pi*(fc)*((0:(N-1))/fs));
    %     case 2; input_signal = bpsk_generator(fs*10^6,fc*10^6,PW*10^3);
    %     case 3; input_signal = qpsk_generator(fs*10^6,fc*10^6,PW*10^3);
    %     case 4; input_signal = lfm_generator(fs) ;
    %     case 5; input_signal = nlfm_generator(fs) ;
    %     otherwise;
    % end
    % x = awgn(input_signal,10);
    % M=128;
    % dalpha=fs/N;
    % df=dalpha*M;
    %-----------Definition of Parameters-----------%
    Np = pow2(nextpow2(fs/df));     %输入信道数
    L = Np/4;                       %抽取因子,同一行连续列的相邻点的偏移,L表示每次滑动的数据点数
    P = pow2(nextpow2(fs/dalpha/L));%信道矩阵的列数,P表示滑动次数
    N = P*L;                        %输入数据的点数
    %----------Input Channelization----------------%
    if length(x) < N
        x(N) = 0;
    elseif length(x) > N
        x = x(1:N);
    end
    NN = (P-1)*L+Np;
    xx = x;
    xx(NN) = 0;
    xx = xx(:);
    X = zeros(Np,P);
    for k = 0:P-1
        X(:,k+1) = xx(k*L+1:k*L+Np);  
        %输入数据xx的1到Np存入X第一列,L+1到L+Np存入第二列,2L+1到2L+Np存入第三列,以此类推。即L为偏移,Np为每段长度。
    end
    %------------------Windowing------------------%
    % a = hamming(Np);
    a=kaiser(Np);
    XW = diag(a)*X;   %每段加窗,减少频率泄露
    % XW = X;         %不加窗,可与加窗的效果作比较
    %---------------第一次傅里叶变换--------------------%
    XF1 = fft(XW);  
    % clear XW;   
    XF1 = fftshift(XF1);   
    XF1 = [XF1(:,P/2+1:P) XF1(:,1:P/2)];  %这两行的功能是把傅里叶变换的结果上半块和下半块交换,左半块和右半块互换
    
    %-------------下变频------------------%
    E = zeros(Np,P);
    
    for k = -Np/2:Np/2-1
        for m = 0:P-1
            E(k+Np/2+1,m+1) = exp(-i*2*pi*k*m*L/Np);
        end
    end
    XD = XF1.*E;   
    % clear XF1; 
    XD = conj(XD');                      %XD转置的复共轭
    % clear ('XF1', 'E', 'XW', 'X', 'x'); 
    %-----------乘积运算--------------------%
    XM = zeros(P,Np^2);
    for k = 1:Np
        for l = 1:Np
            XM(:,(k-1)*Np+l) = (XD(:,k).*conj(XD(:,l)));
        end
    end
    % clear XD;
    %-----------第二次傅里叶变换-----------------------%
    XF2 = fft(XM);  
    % clear XM;
    XF2 = fftshift(XF2);
    XF2 = [XF2(:,Np^2/2+1:Np^2) XF2(:,1:Np^2/2)];
    % length(XF2);
    XF2 = XF2(P/4:3*P/4,:);
    % M = abs(XF2);
     M = XF2;
    % clear XF2;  %Absolute value and complex magnitude
    
    %%%%%%%%%%%%%%%%频率分辨率和循环频率分辨率%%%%%%%%%%
    alphao = (-1:1/N:1)*fs;     %about the variable N!!!
    fo = (-0.5:1/Np:0.5)*fs;
    Sx = zeros(Np+1,2*N+1);     %about the variable N!!!
    for k1 = 1:P/2+1
        for k2 = 1:Np^2
            if rem(k2,Np) == 0
                l = Np/2-1;    
            else
                l = rem(k2,Np)-Np/2-1; 
            end
            k = ceil(k2/Np)-Np/2-1; 
            p = k1-P/4-1;
            alpha = (k-l)/Np+(p-1)/L/P;
            f = (k+l)/2/Np;
            if alpha < -1 || alpha > 1
                k2 = k2+1;
            elseif f < -0.5 || f > 0.5
                k2 = k2+1;
    %         elseif rem(k+l,2)==0 && rem(1+N*(alpha+1),1)==0
            else
                kk = 1+Np*(f+0.5);
                ll = 1+N*(alpha + 1);
                Sx(round(kk), round(ll)) = M(k1,k2); 
            end
        end
    end
    Sx = Sx./max(max(Sx)); % Normalizes the magnitudes of the values in output matrix (maximum = 1)
    

    参考文献

    [1] Gardner W. A . The spectral correlation theory of cyclostationary time-series.
    [2] WILLIAM A. GARDNER -Spectral Correlation of Modulated Signals: Part I-Analog Modulation
    [3] WILLIAM A. GARDNER -Spectral Correlation of Modulated Signals: Part II-Digital Modulation
    文献 [2] [3]给出了深入的介绍性处理,这是充分理解谱相关概念的先决条件。

    展开全文
  • z^\alpha (f) = H(f + \alpha /2){H^*}(f - \alpha /2)S_x^\alpha (f)Szα​(f)=H(f+α/2)H∗(f−α/2)Sxα​(f)(5) 2 周期抽样函数的循环自相关函数和谱相关密度函数 周期抽样函数可以表示为 $g(t) = \sum\limits_{...

    1 引言

    R^xα(f)0\hat R_x^\alpha (f) \equiv 0,对α0,R^x(τ)0\forall \alpha \ne 0,{\hat R_x}(\tau ) \ne 0,则x(t)x(t)为纯平稳信号。

    R^xα(f)0\hat R_x^\alpha (f) \ne 0,当α=mT0\alpha = \frac{m}{{{T_0}}},则x(t)x(t)具备纯循环平稳性,周期为T0{T_0}

    R^xα(f)0\hat R_x^\alpha (f) \ne 0,其中α\alpha不全为1T0\frac{1}{{{T_0}}}的整数倍,则x(t)x(t)为循环平稳的。

    R^yα(τ)=n,m=tr{[R^xα(nm)/T0(τ)ejπ(n+m)τ/T0]rnmα(τ)}\hat R_y^\alpha (\tau ) = \sum\limits_{n,m = - \infty }^\infty {\operatorname{tr} \left\{ {\left[ {\hat R_x^{\alpha - (n - m)/{T_0}}(\tau ){{\text{e}}^{ - j\pi (n + m)\tau /{T_0}}}} \right] \otimes r_{nm}^\alpha ( - \tau )} \right\}}

    rnmα(τ)gn(t+12τ)gm(t12τ)ei2παtdtr_{nm}^\alpha (\tau ) \triangleq \int_{ - \infty }^\infty {{{g'}_n}} \left( {t + \frac{1}{2}\tau } \right)g_m^*\left( {t - \frac{1}{2}\tau } \right){{\text{e}}^{ - {\text{i}}2\pi \alpha t}}{\text{d}}t

    rnmα(τ)S^rα(f)=(gn(t+12τ)gm(t12τ)ej2παtdt)ej2πfτdτr_{nm}^\alpha (\tau ) \leftrightarrow \hat S_r^\alpha (f) = \int_{ - \infty }^\infty {\left( {\int_{ - \infty }^\infty {{{g'}_n}\left( {t + \frac{1}{2}\tau } \right)g_m^*\left( {t - \frac{1}{2}\tau } \right){{\text{e}}^{ - j2\pi \alpha t}}{\text{d}}t} } \right){e^{ - j2\pi f\tau }}d\tau }

    =gn(t+12τ)ej2π(f+α/2)(t+τ/2)d(t+12τ)gm(t12τ)ej2π(fα/2)(tτ/2)d(t12τ)= \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {{{g'}_n}\left( {t + \frac{1}{2}\tau } \right){e^{ - j2\pi (f + \alpha /2)(t + \tau /2)}}{\text{d(}}t + \frac{1}{2}\tau )g_m^*\left( {t - \frac{1}{2}\tau } \right)} {e^{j2\pi (f - \alpha /2)(t - \tau /2)}}d(t - \frac{1}{2}\tau )}

    =gn(t+12τ)ej2π(f+α/2)(t+τ/2)d(t+τ/2)(gm(t12τ)ej2π(fα/2)(tτ/2))d(t12τ)= \int_{ - \infty }^\infty {{{g'}_n}\left( {t + \frac{1}{2}\tau } \right){e^{ - j2\pi (f + \alpha /2)(t + \tau /2)}}{\text{d(}}t + \tau /2)\int_{ - \infty }^\infty {{{\left( {{g_m}\left( {t - \frac{1}{2}\tau } \right){e^{ - j2\pi (f - \alpha /2)(t - \tau /2)}}} \right)}^*}d(t - \frac{1}{2}\tau )} }

    =Gn(f+α/2)Gm(fα/2)= {{G'}_n}(f + \alpha /2)G_m^*(f - \alpha /2)(1)

    S^yα(f) =n,m=Gn(f+12α)S^xα(nm)/T0(fn+m2T0)Gm(f12α)\hat S_y^\alpha (f){\text{ }} = \sum\limits_{n,m = - \infty }^\infty {{G_n}} \left( {f + \frac{1}{2}\alpha } \right)\hat S_x^{\alpha - (n - m)/{T_0}}\left( {f - \frac{{n + m}}{{2{T_0}}}} \right){G'_m}{\left( {f - \frac{1}{2}\alpha } \right)^*}(2)

    信号x(t)x(t)通过线性时不变系统后的输出为

    z(t)=h(t)x(t)=h(u)x(tu)duz(t) = h(t) \otimes x(t) = \int_{ - \infty }^\infty {h(u)x(t - u)du}(3)

    其功率谱为

    Sz(f)=H(f)2Sx(f){S_z}(f) = {\left| {H(f)} \right|^2}{S_x}(f)(4)

    谱相关密度函数为

    Szα(f)=H(f+α/2)H(fα/2)Sxα(f)S_z^\alpha (f) = H(f + \alpha /2){H^*}(f - \alpha /2)S_x^\alpha (f)(5)

    2 周期抽样函数的循环自相关函数和谱相关密度函数

    周期抽样函数可以表示为

    $g(t) = \sum\limits_{n = - \infty }^\infty {\delta (t - nT)} $(6)

    Rxα(τ)limx1TT/2T/2x(t+τ/2)x(tτ/2)dt=x(t+τ/2)x(tτ/2)tR_x^\alpha (\tau ) \triangleq \mathop {\lim }\limits_{x \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {x(t + \tau /2)x(t - \tau /2)dt} = {\left\langle {x(t + \tau /2)x(t - \tau /2)} \right\rangle _t},得

    Rgβ(τ)=n=m=δ(t+τ/2nT)δ(tτ/2mT)ej2πβttR_g^\beta (\tau ) = {\left\langle {\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\delta (t + \tau /2 - nT)\delta (t - \tau /2 - mT){e^{ - j2\pi \beta \cdot t}}} } } \right\rangle _t}

    =n=m=δ((τ/2+mT)+τ/2nT)δ(tτ/2mT)ej2πβ(tτ/2+τ/2)t= {\left\langle {\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\delta ((\tau /2 + mT) + \tau /2 - nT)\delta (t - \tau /2 - mT){e^{ - j2\pi \beta (t - \tau /2 + \tau /2)}}} } } \right\rangle _t}

    =ejπβτn=m=δ(τ+(mn)T)m=δ(tτ/2mT)ej2πβ(tτ/2)t= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\delta (\tau + (m - n)T)} } {\left\langle {\sum\limits_{m = - \infty }^\infty {\delta (t - \tau /2 - mT){e^{ - j2\pi \beta (t - \tau /2)}}} } \right\rangle _t}

    u=tτ/2ejπβτn=m=δ(τ+(mn)T)limT1T(T+τ)/2(Tτ)/2m=δ(umT)ej2πβmTdu{\underline{\underline {u = t - \tau /2}} } {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\delta (\tau + (m - n)T)} } \mathop {\lim }\limits_{T' \to \infty } \frac{1}{{T'}}\int_{ - (T' + \tau )/2}^{(T' - \tau )/2} {\sum\limits_{m = - \infty }^\infty {\delta (u - mT){e^{ - j2\pi \beta mT}}} du}

    =ejπβτn=δ(τnT)limN12NT+2εNTεNT+εm=NNδ(umT)ej2πmmdu= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)} \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2NT + 2\varepsilon }}\int_{ - NT - \varepsilon }^{NT + \varepsilon } {\sum\limits_{m = - N}^N {\delta (u - mT){e^{ - j2\pi m \cdot m}}} du}

    =ejπβτn=δ(τnT)limN12NT+2εNTεNT+εm=NNδ(umT)ej2πmmdu= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)} \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2NT + 2\varepsilon }}\int_{ - NT - \varepsilon }^{NT + \varepsilon } {\sum\limits_{m = - N}^N {\delta (u - mT){e^{ - j2\pi m \cdot m}}} du}

    =ejπβτn=δ(τnT)limN12NT+2εm=NNNTεNT+εδ(umT)du= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)} \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2NT + 2\varepsilon }}\sum\limits_{m = - N}^N {\int_{ - NT - \varepsilon }^{NT + \varepsilon } {\delta (u - mT)du} }

    =ejπβτn=δ(τnT)limN2N+12NT+2ε= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)} \mathop {\lim }\limits_{N \to \infty } \frac{{2N + 1}}{{2NT + 2\varepsilon }}

    =ejπβτn=δ(τnT)= {e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)}(7)

    则抽样函数的谱相关密度函数为

    Sgβ(f)=Rgβ(τ)ej2πfτdτS_g^\beta (f) = \int_{ - \infty }^\infty {R_g^\beta (\tau ){e^{ - j2\pi f\tau }}d\tau }

    =ejπβτn=δ(τnT)ej2πfτdτ= \int_{ - \infty }^\infty {{e^{ - j\pi \beta \tau }}\sum\limits_{n = - \infty }^\infty {\delta (\tau - nT)} {e^{ - j2\pi f\tau }}d\tau }

    =n=ejπβnTej2πfnTδ(τnT)dτ= \sum\limits_{n = - \infty }^\infty {{e^{ - j\pi \beta nT}}{e^{ - j2\pi fnT}}\int_{ - \infty }^\infty {\delta (\tau - nT)d\tau } }

    =n=ej2πnT(f+β2)= \sum\limits_{n = - \infty }^\infty {{e^{ - j2\pi nT(f + \frac{\beta }{2})}}}

    =n=ej2πTn(f+β2)= \sum\limits_{n = - \infty }^\infty {{e^{ - j\frac{{2\pi }}{{T'}}n(f + \frac{\beta }{2})}}}(8)

    其中T=1/TT' = 1/T,由n=δ(tnT)=1Tn=ej2πTnt\sum\limits_{n = - \infty }^\infty {\delta (t - nT)} = \frac{1}{T}\sum\limits_{n = - \infty }^\infty {{e^{j\frac{{2\pi }}{T}n \cdot t}}},得

    Sgβ(f)=Tk=δ(f+β2kT)=1Tk=δ(f+β2kT)S_g^\beta (f) = T'\sum\limits_{k = - \infty }^\infty {\delta (f + \frac{\beta }{2} - kT')}= \frac{1}{T}\sum\limits_{k = - \infty }^\infty {\delta (f + \frac{\beta }{2} - \frac{k}{T})}(9)

    3 信号相乘前后谱相关密度函数关系

    假设x(t)=r(t)s(t)x(t) = r(t)s(t),由

    Rx(t,τ)=E[x(t+τ2)x(tτ2)]=Rr(t,τ)Rs(t,τ){R_x}(t,\tau ) = E[x(t + \frac{\tau }{2}){x^*}(t - \frac{\tau }{2})] = {R_r}(t,\tau ){R_s}(t,\tau )(10)

    ${R_x}(t,\tau ) = \sum\limits_\alpha {R_x^\alpha (\tau ){e^{j2\pi \alpha t}}} $(11)

    x(t)x(t)的循环自相关函数为

    Rxα(τ)=limT1TT/2T/2Rx(t,τ)ej2παtdtR_x^\alpha (\tau ) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {{R_x}(t,\tau ){e^{ - j2\pi \alpha t}}dt}

    =limT1TT/2T/2Rr(t,τ)Rs(t,τ)ej2παtdt= \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {{R_r}(t,\tau ){R_s}(t,\tau ){e^{ - j2\pi \alpha t}}dt}

    =limT1TT/2T/2βRr(t,τ)ej2π(αβ)tRs(t,τ)ej2πβtdt= \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {\sum\limits_\beta {{R_r}(t,\tau ){e^{ - j2\pi (\alpha - \beta )t}}{R_s}(t,\tau ){e^{ - j2\pi \beta t}}} dt}

    =βlimT1TT/2T/2Rr(t,τ)ej2π(αβ)tRs(t,τ)ej2πβtdt= \sum\limits_\beta {\mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {{R_r}(t,\tau ){e^{ - j2\pi (\alpha - \beta )t}}{R_s}(t,\tau ){e^{ - j2\pi \beta t}}dt} }

    =βRrαβ(τ)Rsβ(τ)= \sum\limits_\beta {R_r^{\alpha - \beta }(\tau )R_s^\beta (\tau )}

    =βRrβ(τ)Rsαβ(τ)= \sum\limits_\beta {R_r^\beta (\tau )R_s^{\alpha - \beta }(\tau )} (12式)

    由(12)可知,两个信号相乘后的循环自相关函数等于两个循环自相关函数在离散$\alpha $域的卷积。

    其谱相关密度函数为

    Sxα(f)=Rxα(τ)ej2πfτdτS_x^\alpha (f) = \int_{ - \infty }^\infty {R_x^\alpha (\tau ){e^{ - j2\pi f\tau }}d\tau }

    =βRrβ(τ)Rsαβ(τ)ej2πfτdτ= \int_{ - \infty }^\infty {\sum\limits_\beta {R_r^\beta (\tau )R_s^{\alpha - \beta }(\tau )} {e^{ - j2\pi f\tau }}d\tau }

    =βRrβ(τ)ej2πFτRsαβ(τ)ej2π(fF)τdτdF= \sum\limits_\beta {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {R_r^\beta (\tau ){e^{ - j2\pi F\tau }}R_s^{\alpha - \beta }(\tau ){e^{ - j2\pi (f - F)\tau }}d\tau dF} } }

    =βRrβ(τ)ej2πFτRsαβ(τ)ej2π(fF)τdτdF= \sum\limits_\beta {\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {R_r^\beta (\tau ){e^{ - j2\pi F\tau }}R_s^{\alpha - \beta }(\tau ){e^{ - j2\pi (f - F)\tau }}d\tau dF} } }

    =βSrβ(F)Ssαβ(fF)dF= \int_{ - \infty }^\infty {\sum\limits_\beta {S_r^\beta (F)S_s^{\alpha - \beta }(f - F)} dF}

    =βSrβ(f)Ssαβ(f)= \sum\limits_\beta {S_r^\beta (f) \otimes S_s^{\alpha - \beta }(f)} (13式)

    其中\otimes表示卷积,由此可知两个信号相乘的谱相关等于两个谱相关在连续ff域和离散α\alpha域的双重卷积。

    4 周期信号(Almost periodic)的循环自相关与谱相关密度函数

    p(t)p(t)为周期信号,其傅里叶级数可以表示为

    p(t)=βPβej2πβtp(t) = \sum\limits_\beta {{P_\beta }{e^{j2\pi \beta t}}} (14式)

    其中β=nT\beta = \frac{n}{T}TT为周期,则

    Pβ=p(t)ej2πβtt{P_\beta } = {\left\langle {p(t){e^{ - j2\pi \beta t}}} \right\rangle _t} (15式)

    由循环自相关函数的定义得

    Rpα(τ)=p(t)p(tτ)ej2πα(tτ/2)tR_p^\alpha (\tau ) = {\left\langle {p(t){p^*}(t - \tau ){e^{ - j2\pi \alpha (t - \tau /2)}}} \right\rangle _t}

    =βp(t)ej2πβtp(tτ)ej2π(αβ)(tτ)ejπ(2βα)τt= {\left\langle {\sum\limits_\beta {p(t){e^{ - j2\pi \beta t}}{p^*}(t - \tau ){e^{ - j2\pi (\alpha - \beta )(t - \tau )}} \cdot {e^{j\pi (2\beta - \alpha )\tau }}} } \right\rangle _t}

    =βPβPαβejπ(2βα)τ= \sum\limits_\beta {{P_\beta }P_{\alpha - \beta }^*{e^{j\pi (2\beta - \alpha )\tau }}} (16式)

    则其谱相关密度函数为

    Spα(f)=βPβPαβejπ(2βα)τej2πfτdτS_p^\alpha (f) = \int_{ - \infty }^\infty {\sum\limits_\beta {{P_\beta }P_{\alpha - \beta }^*{e^{j\pi (2\beta - \alpha )\tau }}} {e^{ - j2\pi f\tau }}d\tau }

    =βPβPαβej2π(fβ+α/2)τdτ= \sum\limits_\beta {{P_\beta }P_{\alpha - \beta }^*\int_{ - \infty }^\infty {{e^{ - j2\pi (f - \beta + \alpha /2)\tau }}d\tau } }

    =βPβPαβδ(fβ+α/2)= \sum\limits_\beta {{P_\beta }P_{\alpha - \beta }^*\delta (f - \beta + \alpha /2)} (17式)

    5 连续平稳信号被理想抽样后的谱相关密度函数

    连续平稳信号x(t)x(t)被理想抽样后的信号可以表示为

    y(t)=x(t)g(t)=n=x(nT)δ(tnT)y(t) = x(t)g(t) = \sum\limits_{n = - \infty }^\infty {x(nT)\delta (t - nT)} (18式)

    其中,g(t)=n=δ(tnT)g(t) = \sum\limits_{n = - \infty }^\infty {\delta (t - nT)},则由(13)得y(t)y(t)的谱相关密度函数为

    Syα(f)=βSgβ(F)Sxαβ(fF)dFS_y^\alpha (f) = \int_{ - \infty }^\infty {\sum\limits_\beta {S_g^\beta (F)S_x^{\alpha - \beta }(f - F)} dF}

    =β1T2n=δ(F+β2nT)Sxαβ(fF)dF= \sum\limits_\beta {\int_{ - \infty }^\infty {\frac{1}{{{T^2}}}\sum\limits_{n = - \infty }^\infty {\delta (F + \frac{\beta }{2} - \frac{n}{T})} S_x^{\alpha - \beta }(f - F)dF} }

    =m=1T2n=δ(F+β2nT)Sxαβ(fF)dF= \sum\limits_{m = - \infty }^\infty {\int_{ - \infty }^\infty {\frac{1}{{{T^2}}}\sum\limits_{n = - \infty }^\infty {\delta (F + \frac{\beta }{2} - \frac{n}{T})} S_x^{\alpha - \beta }(f - F)dF} }

    =1T2m=n=SxαmT(f+m2TnT)δ(F+m2TnT)dF= \frac{1}{{{T^2}}}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {S_x^{\alpha - \frac{m}{T}}(f + \frac{m}{{2T}} - \frac{n}{T})} \int_{ - \infty }^\infty {\delta (F + \frac{m}{{2T}} - \frac{n}{T})dF} }

    =1T2m=n=SxαmT(f+m2TnT)= \frac{1}{{{T^2}}}\sum\limits_{m = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {S_x^{\alpha - \frac{m}{T}}(f + \frac{m}{{2T}} - \frac{n}{T})} } (19式)

    6 连续平稳信号与其周期采样序列的谱相关密度函数关系

    连续平稳信号x(t)非对称形式的循环自相关函数为
    在这里插入图片描述(1)
    在这里插入图片描述
    由(1)定义离散时间序列x(nT0)x(nT_0)的循环自相关函数为
    在这里插入图片描述(2)
    x(t)的谱相关密度函数为
    在这里插入图片描述(3)
    由(3)定义x(nT0)x(nT_0)的谱相关密度函数为
    在这里插入图片描述(4)
    将恒等式
    在这里插入图片描述(5)
    代入(2),得
    R~xα(kT0)limN12N+1n=NNx(nT0+kT0)x(nT0)ej2πα(n+k/2)T0\tilde R_x^\alpha (k{T_0}) \triangleq \mathop {\lim }\limits_{N \to \infty } \frac{1}{{2N + 1}}\sum\limits_{n = - N}^N {x(n{T_0} + k{T_0})x(n{T_0})} {e^{ - j2\pi \alpha (n + k/2){T_0}}}

    nT0um=limT1TT/2T/2x(u+kT0)x(u)ej2πα(u+kT0/2)ej2πmu/T0du{\underline{\underline {n{T_0} \to u}} } \sum\limits_{m = - \infty }^\infty {\mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {x(u + k{T_0})x(u)} } {e^{ - j2\pi \alpha (u + k{T_0}/2)}}{e^{ - j2\pi mu/{T_0}}}du

    =m=limT1TT/2T/2x(u+kT0)x(u)ej2πα(u+kT0/2)ej2πm/T0(u+kT0/2)ejπmkdu= \sum\limits_{m = - \infty }^\infty {\mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {x(u + k{T_0})x(u)} } {e^{ - j2\pi \alpha (u + k{T_0}/2)}}{e^{ - j2\pi m/{T_0}(u + k{T_0}/2)}}{e^{j\pi mk}}du

    =m=limT1TT/2T/2x(u+kT0)x(u)ej2π(α+m/T0)(u+kT0/2)duejπmk= \sum\limits_{m = - \infty }^\infty {\mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} {x(u + k{T_0})x(u)} } {e^{ - j2\pi (\alpha + m/{T_0})(u + k{T_0}/2)}}du \cdot {e^{j\pi mk}}

    =m=R^xα + m/T0(kT0)ejπmk= \sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{j\pi mk}} (6式)

    将(6)代入(4),得一种错误的结果如下

    S~xα(f)=k=[m=R^xα + m/T0(kT0)ejπmk]ej2πkT0f\tilde S_x^\alpha (f) = \sum\limits_{k = - \infty }^\infty {[\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{j\pi mk}}]{e^{ - j2\pi k{T_0}f}}}

    =k=m=R^xα + m/T0(kT0)ej2πkT0(fm2T0)= \sum\limits_{k = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{ - j2\pi k{T_0}(f - \frac{m}{{2{T_0}}})}}}

    kT0τ1T0m=R^xα + m/T0(τ)ej2π(fm2T0)τdτ{\underline{\underline {k{T_0} \to \tau }} } \frac{1}{{{T_0}}}\int_{ - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(\tau ){e^{ - j2\pi (f - \frac{m}{{2{T_0}}})\tau }}} d\tau }

    =1T0m=R^xα + m/T0(τ)ej2π(fm2T0)τdτ= \frac{1}{{{T_0}}}\sum\limits_{m = - \infty }^\infty {\int_{ - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(\tau ){e^{ - j2\pi (f - \frac{m}{{2{T_0}}})\tau }}d\tau } }

    =1T0m=S^xα + m/T0(fm2T0)= \frac{1}{{{T_0}}}\sum\limits_{m = - \infty }^\infty {\hat S_x^{\alpha {\text{ + }}m/{T_0}}(f - \frac{m}{{2{T_0}}})} (20式)

    正确结果应为Sxα(f)S_x^\alpha (f)

    S~xα(f)=k=[m=R^xα + m/T0(kT0)ejπmk]ej2πkT0f\tilde S_x^\alpha (f) = \sum\limits_{k = - \infty }^\infty {[\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{j\pi mk}}]{e^{ - j2\pi k{T_0}f}}}

    =k=m=R^xα + m/T0(kT0)ej2πkT0(fm2T0)= \sum\limits_{k = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{ - j2\pi k{T_0}(f - \frac{m}{{2{T_0}}})}}}

    =k=n=m=R^xα + m/T0(kT0)ej2πkT0(fm2T0)ej2πnk= \sum\limits_{k = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{ - j2\pi k{T_0}(f - \frac{m}{{2{T_0}}})}} \cdot {e^{j2\pi nk}}} }

    =k=n=m=R^xα + m/T0(kT0)ej2πkT0(fm2T0nT0)= \sum\limits_{k = - \infty }^\infty {\sum\limits_{n = - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(k{T_0})} {e^{ - j2\pi k{T_0}(f - \frac{m}{{2{T_0}}} - \frac{n}{{{T_0}}})}}} }

    kT0τ1T0m=R^xα + m/T0(τ)ej2π(fm2T0nT0)τdτ{\underline{\underline {k{T_0} \to \tau }} } \frac{1}{{{T_0}}}\int_{ - \infty }^\infty {\sum\limits_{m = - \infty }^\infty {\hat R_x^{\alpha {\text{ + }}m/{T_0}}(\tau ){e^{ - j2\pi (f - \frac{m}{{2{T_0}}} - \frac{n}{{{T_0}}})\tau }}} d\tau }