小波变换估计载波速率整理.docx

1、小波变换估计载波速率整理小波变换估计载波速率作者： 时间：2011-9-16已知：接收机采集的数字调制信号，包括MASK、MFSK、MPSK、MQAM估计：码速率 （估计的是符号速率，码片速率=符号速率码元宽度，码元宽度为一个符号的二进制位宽，如8PSK，8=23，所以码元宽度为3）参考文章：1、2、-博士论文3、matlab中cwt函数源码及武瑞娟-Haar小波详解及matlab源码.doc4、一种用于信号侦测的码元速率估计方法论文1、方法概述采用小波变换+自相关 估计数字调制码速率；2、原理概述单极性脉冲序列的功率谱在码速率的整数倍处存在离散谱线，检测这些谱线即可实现码速率估计（【2】的2

2、.1节）。因此对于数字调制信号的基带信号构建单极性波形后即可估计码速率。数字基带信号波形在码元变换处存在突跳，小波变换能有效检测信号的奇异点，同时小波变换的模值序列可构成单极性随机脉冲波形。当待分析的信号与小波基正确匹配时，信号的主要能量都分布在少数系数上。若把这些系数的模值看成是单极性脉冲序列，则其功率谱包含了信号码速率信息。参见参考文档2的2.3节。构造单极性脉冲序列：1）对于ASK，即幅度调制的信号，求其包络，然后求包络的差分，在码元变换出存在突跳，构造除了单极性脉冲序列；2）对于FSK，码元变化表现为瞬时频率变化，因此求其瞬时频率fn，然后瞬时频率的差分构造单极性脉冲；3）对于PSK，

3、码元变化表现为瞬时相位的变化，因此求其瞬时相位ph，然后瞬时相位的差分构造单极性脉冲；4）对于QAM，码元变化即表现为幅度变化，又表现为相位变化，选择瞬时相位ph的差分构造单极性脉冲；瞬时参数的求法：1）瞬时幅度 2）瞬时相位 ，修正为无折叠相位，C（0）=0；3) 瞬时频率 3、码速率估计步骤Fdcompute=2*fs/(upindex-lowindex)；其中fs为采样率，upindex和lowindex分别为主峰两侧峰值的索引如下图所示 其中，haar小波变换未使用matlab自带的cwt，而是等效函数，可以方便的改写为c+等其它语言。附件1、码速率估计matlab仿真%psk qam

4、 fsk ask可以测出码速率clc;clear;j=32;fc=1000;fd=120;fs=1200;select=menu(调制方式,2PSK,4PSK,8PSK,2QAM,4QAM,8QAM, 2FSK,4FSK, 2ASK,4ASK);switch selectcase 1, M=2; x=0 1 1 0 0 1 0 1 0 1 0 1 1 0 1 0 0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,psk,M);%psk qam fsk pam case 2, M=4; x=0 2 1 3 0

5、1 2 3 2 1 0 1 1 2 1 0 2 3 2 1 3 2 0 1 3 0 2 0 1 3 2 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,psk,M);%psk qam fsk pam case 3, M=8; x=0 1 5 3 7 5 6 0 3 5 0 4 3 2 6 0 6 3 4 5 7 2 0 1 3 1 6 4 1 3 2 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,psk,M);%psk qam fsk pam case 4, M=2; x=0 1 1 0 0 1 0 1 0 1 0 1 1 0 1 0

6、0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,qam,M);%psk qam fsk pam case 5, M=4; x=0 2 1 2 3 2 1 3 2 1 0 1 1 2 1 0 2 3 2 1 3 2 0 1 2 3 2 1 1 3 2 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,qam,M);%psk qam fsk pam case 6, M=8; x=0 1 2 3 4 5 6 7 3 2 0 4 3 3 6 0 0 3 0 5 7 2 0 1 3 4 6

7、5 1 3 2 0;%randint(1,j); sn,t=dmod(x,fc,fd,fs,qam,M);%psk qam fsk pam case 7, %2fsk M=2; x=0 1 0 1 1 0 0 1 0 1 0 1 0 1 1 0 1 0 0 1 0 1 0 1 0 0 1 0 1 0 1 0;%randint(1,j); ts=1/fs; N=fix(fs/fd); t=0:ts:j*N*ts; m1=1*cos(2*pi*fc*t);%wc=2*pi*fc 是cos 是pi的整数倍，因此取值只有0和1 m2=1*cos(2*pi*2*fc*t); %恒为1 for i=1:j

8、 if x(i)=1; for k=1:N; sn(i-1)*N+k)=x(i)*m1(k); end; elseif x(i)=0; for k=1:N; sn(i-1)*N+k)=(1-x(i)*m2(k); end; end; end;% sn,t=dmod(x,fc,fd,fs,fsk,M);%psk qam fsk pam 连续相位不可以测case 8, %4fsk M=4; x=0 2 1 3 0 2 0 1 3 2 0 1 1 3 1 0 0 2 0 1 3 2 0 1 3 0 2 0 1 3 2 0;%randint(1,j); ts=1/fs; N=fix(fs/fd); t

9、=0:ts:j*N*ts; m1=cos(2*pi*fc*t); m2=cos(2*pi*1.5*fc*t); m3=cos(2*pi*1.9*fc*t); m4=cos(2*pi*2.5*fc*t); for i=1:j if x(i)=0; for k=1:N sn(i-1)*N+k)=(1-x(i)*m1(k); end elseif x(i)=1; for k=1:N sn(i-1)*N+k)=x(i)*m2(k); end elseif x(i)=2; for k=1:N sn(i-1)*N+k)=(x(i)-1)*m3(k); end elseif x(i)=3; for k=1:

10、N sn(i-1)*N+k)=(x(i)-2)*m4(k); end end end;% M=4;% x=0 2 1 3 0 2 0 1 3 2 0 1 1 3 1 0 0 2 0 1 3 2 0 1 3 0 2 0 1 3 2 0;%randint(1,j);% sn,t=dmod(x,fc,fd,fs,fsk,M);%psk qam fsk pam case 9, %2Ask M=2; x=0 1 1 0 0 1 0 1 0 1 0 1 1 0 1 0 0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0;%randint(1,j); ts=1/fs; N=fix(fs/fd);

11、 t=0:ts:j*N*ts; m=1*cos(2*pi*fc*t);%wc=2*pi*fc 是cos 是pi的整数倍，因此取值只有0和1 for i=1:j for k=1:N sn(i-1)*N+k)=x(i)*m(k); %2ASK信号可以表示为一个单极性脉冲与一个正弦载波相乘 y=1*20480 即M*N个 end end % sn,t=dmod(x,fc,fd,fs,ask,M);%ask case 10, %4Ask M=4; x=0 2 1 3 0 2 0 1 3 2 0 1 1 3 1 0 0 2 0 1 3 2 0 1 3 0 2 0 1 3 2 0;%randint(1,j

12、); ts=1/fs; N=fix(fs/fd); t=0:ts:j*N*ts; m=1*cos(2*pi*fc*t); for i=1:j if x(i)=0; for k=1:N sn(i-1)*N+k)=x(i)*m(k); end elseif x(i)=1; for k=1:N sn(i-1)*N+k)=x(i)*m(k); end elseif x(i)=2; for k=1:N sn(i-1)*N+k)=x(i)*m(k); end elseif x(i)=3; for k=1:N sn(i-1)*N+k)=x(i)*m(k); end end endend; %M=4;s=aw

13、gn(sn,15,measured,db);%snr9 t=1/fs:1/fs:j/fd;DataSrc=menu(数据源,生成,读取);switch DataSrccase 2, fid = fopen(E:autoDemfsk200886 AM 103805 M 166.775,50IQ.dat,r); index=1*2*8192; Len = 8192; s = fread(fid, Len*2+14+index, short); fclose(fid); fs=1000000;%采样率 s=s(15+index:Len*2+14+index); subplot(4,1,1);plot

14、(s); II=s(1:2:Len*2); QQ=s(2:2:Len*2); subplot(4,1,2);plot(II); subplot(4,1,3);plot(QQ); s=II;case 1, subplot(4,1,1);plot(s); tmps=hilbert(s); II=real(tmps); QQ=imag(tmps); subplot(4,1,2);plot(II); subplot(4,1,3);plot(QQ); %s=II;end;N=length(s);% if (select = 1 | select = 2 | select = 3 | select = 4

15、 | select = 5 | select = 6| select=7 | select = 8); %psk FSKselect = 1 | select = 2 | select = 3 | %求相位 for i=1:N; phase(i)= atan2(QQ(i), II(i); end; pi=3.1415926; %修正序列 tmpPhaseCorr(1)=0; for i =1:N-1; if (phase(i+1)-phase(i) pi) tmpPhaseCorr(i+1) = tmpPhaseCorr(i) - 2 * pi; elseif(phase(i)-phase(i

16、+1) pi) tmpPhaseCorr(i+1) = tmpPhaseCorr(i) + 2 * pi; else tmpPhaseCorr(i+1) = tmpPhaseCorr(i); end; end; phase = phase+tmpPhaseCorr; phase = FIT(phase); fn = phase(2:end) - phase(1:end-1); for i=1:length(fn); while(fn(i) pi) fn(i)=-(2*pi-fn(i); end; while(fn(i) -pi) fn(i)=2*pi+fn(i); end;% fn(i) =

17、fn(i) * fs/(2*pi); end; if (select = 1 | select = 2 | select = 3 | select = 4 | select = 5 | select = 6) s=fn; %相位差 else %平滑 window = 3; fn = windowsmooth(fn, window); s=(fn(2:end)-fn(1:end-1); %频率差 end; subplot(4,1,3);plot(fn); subplot(4,1,4);plot(s); N=length(s); end;if ( select = 9 | select = 10)

18、; %ASKselect = 4 | for i=1:N; am(i)=sqrt(II(i)*II(i)+QQ(i)*QQ(i); end; %平滑 window = 3; am = windowsmooth(am, window); s=am(2:end)-am(1:end-1); subplot(4,1,3);plot(am); subplot(4,1,4);plot(s); N=length(s);end;ppz=10*log10(fftshift(abs(fft(sn);% subplot(2,1,2);plot(ppz);if (select = 1 | select = 2 | s

19、elect = 3); NNN=20;else NNN = 50;end;for i=1:NNN; scale(i)=1+0.2*i;%设置尺度end; s_c_h1 = myHaarcwt(s, scale);%自己写的haar连续小波变换，参照cwt cc=abs(s_c_h1(1,:); for i=2:NNN; cc=cc+abs(s_c_h1(i,:); end;figure(2); %以Mexican hat小波为例%归一化maxv=max(cc);cc=cc/maxv;subplot(411);plot(cc);title(尺度为1.4);%连续小波xx =cc.*cc;subp

20、lot(413);plot(xx);title(2次方);yy = xx.*xx;subplot(414);plot(yy);title(四次方);figure(3);r = xcorr(cc);subplot(3,1,1); plot(abs(r);title(相关结果);rr = xcorr(xx);subplot(3,1,2); plot(abs(rr);title(2次方相关);rrrr = xcorr(yy);subplot(3,1,3); plot(abs(rrrr);title(四次方相关);tmpsignal=r;figure(4)% window = 3;%tmpsignal

21、 = windowsmooth(tmpsignal, window);%subplot(4,1,1); plot(tmpsignal);title(平滑的相关结果);tmpsignal = FindSignal(tmpsignal);subplot(4,1,2); plot(tmpsignal);title(降噪后的相关结果);Find1 = Smoothness(tmpsignal);subplot(4,1,3); plot(Find1);title(1查找信号);Find2 = Smoothness(Find1);subplot(4,1,4); plot(Find2);title(2查找信

22、号);ls=length(Find1);maxindex=-1;max=0;for i=1:ls if (maxFind1(i) max=Find1(i); maxindex=i; end;end;tmprr=Find1;bflag=0;startflag=0;len=0;for i=maxindex+1:ls if tmprr(i)0 if (startflag=0) upindex=i; len=len+1; startflag=1; else len=len+1; end; end;end; upindex=floor(upindex+len/2);%向下取整bflag=0;startf

23、lag=0;len=0;for i=maxindex-1:-1:1; if tmprr(i)0; if (startflag=0) lowindex=i; len=len+1; startflag=1; else len=len+1; end; end;end;lowindex=round(lowindex-len/2);%向上取整tmpfd=2*fs/(upindex-lowindex) function signal = windowsmooth(rudedata, windownum) Ns = length(rudedata); for i =1:Ns; tmpavg=0; if (i

24、 = windownum) for j=i:-1:i-windownum+1; tmpavg = tmpavg+rudedata(j); end; tmpavg=tmpavg/windownum; else for j=i:-1:1; tmpavg = tmpavg+rudedata(j); end; tmpavg=tmpavg/i; end rudedata(i)=tmpavg; end; signal = rudedata;%最小二乘法去线性%procedure FIT(var DATA: array of Double; Ns: integer);function DATA = FIT(

25、sig) SX = 0; SY = 0; Ns = length(sig); for i = 1:Ns; SX = SX+i; SY = SY+sig(i); end; SXOSS = SX/Ns; ST2 = 0; B = 0; for i=1:Ns; T=i-SXOSS; ST2=ST2+T*T; B=B+T*sig(i); end; B = B/ST2; A = (SY - SX*B)/Ns; for i =1:Ns DATA(i)=sig(i)-A-B*i; end;%去除背噪函数function signal = FindSignal(rudedata) Ns = length(ru

26、dedata); maxv=max(rudedata); k = 0; for i = 2:Ns-1; if (rudedata(i) rudedata(i-1) & rudedata(i) rudedata(i+1) peakindex(k + 1) = i; peak(k+1)=rudedata(i); k = k+1; end; end; for i = 2:k; minv=peak(i); for j = peakindex(i-1)+1:peakindex(i)-1; if (rudedata(j) = minv) minv = rudedata(j); valley(i-1) =

27、minv; valleyindex(i-1) = j; end; end; end; for i = 2:k-1; maxv = max(valley(i-1),valley(i); for j = valleyindex(i-1):valleyindex(i); rudedata(j) = rudedata(j) - maxv; if rudedata(j) rudedata(i - 1) & rudedata(i) rudedata(i + 1); peak(k+1)=rudedata(i); sumv = sumv + rudedata(i); k = k + 1; end; end; if (k 0) avg = sumv/k; for i=1:Ns; if (abs(rudedata(i) = avg); signal(i)=0; else signal(i) = rudedata(i); end; end; end;附件2、连续Haar小波变换源码function wcoefs = myHaarcwt(