小波变换估计载波速率整理.docx
《小波变换估计载波速率整理.docx》由会员分享，可在线阅读，更多相关《小波变换估计载波速率整理.docx（20页珍藏版）》请在冰豆网上搜索。
小波变换估计载波速率整理
小波变换估计载波速率
作者：
时间：
2011-9-16
已知：
接收机采集的数字调制信号，包括MASK、MFSK、MPSK、MQAM
估计：
码速率（估计的是符号速率，码片速率=符号速率×码元宽度，码元宽度为一个符号的二进制位宽，如8PSK，8=2^3，所以码元宽度为3）
参考文章：
1、<基于小波变换的移相键控信号符号速率估计（matlab仿真）_西电大宝cc_新浪博客.htm>
2、<数字通信信号的自动识别与参数估计研究（博士）.kdh>--博士论文
3、matlab中cwt函数源码及《武瑞娟-Haar小波详解及matlab源码.doc》
4、《一种用于信号侦测的码元速率估计方法》论文
1、方法概述
采用小波变换+自相关估计数字调制码速率；
2、原理概述
单极性脉冲序列的功率谱在码速率的整数倍处存在离散谱线，检测这些谱线即可实现码速率估计（【2】的2.1节）。
因此对于数字调制信号的基带信号构建单极性波形后即可估计码速率。
数字基带信号波形在码元变换处存在突跳，小波变换能有效检测信号的奇异点，同时小波变换的模值序列可构成单极性随机脉冲波形。
当待分析的信号与小波基正确匹配时，信号的主要能量都分布在少数系数上。
若把这些系数的模值看成是单极性脉冲序列，则其功率谱包含了信号码速率信息。
参见参考文档2的2.3节。
构造单极性脉冲序列：
1）对于ASK，即幅度调制的信号，求其包络
，然后求包络的差分
，在码元变换出存在突跳，构造除了单极性脉冲序列；
2）对于FSK，码元变化表现为瞬时频率变化，因此求其瞬时频率fn，然后瞬时频率的差分构造单极性脉冲；
3）对于PSK，码元变化表现为瞬时相位的变化，因此求其瞬时相位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仿真
%pskqamfskask可以测出码速率
clc;
clear;
j=32;
fc=1000;
fd=120;
fs=1200;
select=menu（'调制方式','2PSK','4PSK','8PSK','2QAM','4QAM','8QAM','2FSK','4FSK','2ASK','4ASK'）;
switchselect
case1,
M=2;
x=[01100101010110100101010100101010];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'psk',M）;%pskqamfskpam
case2,
M=4;
x=[02130123210112102321320130201320];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'psk',M）;%pskqamfskpam
case3,
M=8;
x=[01537560350432606345720131641320];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'psk',M）;%pskqamfskpam
case4,
M=2;
x=[01100101010110100101010100101010];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'qam',M）;%pskqamfskpam
case5,
M=4;
x=[02123213210112102321320123211320];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'qam',M）;%pskqamfskpam
case6,
M=8;
x=[01234567320433600305720134651320];%%randint（1,j）;
[sn,t]=dmod（x,fc,fd,fs,'qam',M）;%pskqamfskpam
case7,%2fsk
M=2;
x=[01011001010101101001010100101010];%%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
fori=1:
j
ifx（i）==1;
fork=1:
N;
sn（（i-1）*N+k）=x（i）*m1（k）;
end;
elseifx（i）==0;
fork=1:
N;
sn（（i-1）*N+k）=（1-x（i））*m2（k）;
end;
end;
end;
%[sn,t]=dmod（x,fc,fd,fs,'fsk',M）;%pskqamfskpam连续相位不可以测
case8,%4fsk
M=4;
x=[02130201320113100201320130201320];%%randint（1,j）;
ts=1/fs;
N=fix（fs/fd）;
t=[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）;
fori=1:
j
ifx（i）==0;
fork=1:
N
sn（（i-1）*N+k）=（1-x（i））*m1（k）;
end
elseifx（i）==1;
fork=1:
N
sn（（i-1）*N+k）=x（i）*m2（k）;
end
elseifx（i）==2;
fork=1:
N
sn（（i-1）*N+k）=（x（i）-1）*m3（k）;
end
elseifx（i）==3;
fork=1:
N
sn（（i-1）*N+k）=（x（i）-2）*m4（k）;
end
end
end;
%M=4;
%x=[02130201320113100201320130201320];%%randint（1,j）;
%[sn,t]=dmod（x,fc,fd,fs,'fsk',M）;%pskqamfskpam
case9,%2Ask
M=2;
x=[01100101010110100101010100101010];%%randint（1,j）;
ts=1/fs;
N=fix（fs/fd）;
t=[0:
ts:
j*N*ts];
m=1*cos（2*pi*fc*t）;%wc=2*pi*fc是cos是pi的整数倍，因此取值只有0和1
fori=1:
j
fork=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
case10,%4Ask
M=4;
x=[02130201320113100201320130201320];%%randint（1,j）;
ts=1/fs;
N=fix（fs/fd）;
t=[0:
ts:
j*N*ts];
m=1*cos（2*pi*fc*t）;
fori=1:
j
ifx（i）==0;
fork=1:
N
sn（（i-1）*N+k）=x（i）*m（k）;
end
elseifx（i）==1;
fork=1:
N
sn（（i-1）*N+k）=x（i）*m（k）;
end
elseifx（i）==2;
fork=1:
N
sn（（i-1）*N+k）=x（i）*m（k）;
end
elseifx（i）==3;
fork=1:
N
sn（（i-1）*N+k）=x（i）*m（k）;
end
end
end
end;
%M=4;
s=awgn（sn,15,'measured','db'）;%snr>9
t=1/fs:
1/fs:
j/fd;
DataSrc=menu（'数据源','生成','读取'）;
switchDataSrc
case2,
fid=fopen（'E:
\autoDem\fsk\200886AM103805M166.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（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;
case1,
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||select==5||select==6||select==7||select==8）;%pskFSKselect==1||select==2||select==3||
%求相位
fori=1:
N;
phase（i）=atan2（QQ（i）,II（i））;
end;
pi=3.1415926;
%修正序列
tmpPhaseCorr
（1）=0;
fori=1:
N-1;
if（phase（i+1）-phase（i）>pi）
tmpPhaseCorr（i+1）=tmpPhaseCorr（i）-2*pi;
elseif（phase（i）-phase（i+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）;
fori=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）=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）;%ASKselect==4||
fori=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||select==3）;
NNN=20;
else
NNN=50;
end;
fori=1:
NNN;
scale（i）=1+0.2*i;%设置尺度
end;
s_c_h1=myHaarcwt（s,scale）;%自己写的haar连续小波变换，参照cwt
cc=abs（s_c_h1（1,:
））;
fori=2:
NNN;
cc=cc+abs（s_c_h1（i,:
））;
end;
figure
（2）;%以Mexicanhat小波为例
%归一化
maxv=max（cc）;
cc=cc/maxv;
subplot（411）;plot（cc）;title（'尺度为1.4'）;%连续小波
xx=cc.*cc;
subplot（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=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查找信号'）;
ls=length（Find1）;
maxindex=-1;
max=0;
fori=1:
ls
if（maxmax=Find1（i）;
maxindex=i;
end;
end;
tmprr=Find1;
bflag=0;
startflag=0;
len=0;
fori=maxindex+1:
ls
iftmprr（i）<=0
ifbflag==0
bflag=1;%startfindupmax
end;
ifstartflag==1
break;
end;
end;
if（bflag==0）
continue;
end;
iftmprr（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;
startflag=0;
len=0;
fori=maxindex-1:
-1:
1;
iftmprr（i）<=0;
ifbflag==0;
bflag=1;%startfindupmax
end;
ifstartflag==1;
break;
end;
end;
if（bflag==0）;
continue;
end;
iftmprr（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）;
fori=1:
Ns;
tmpavg=0;
if（i>=windownum）
forj=i:
-1:
i-windownum+1;
tmpavg=tmpavg+rudedata（j）;
end;
tmpavg=tmpavg/windownum;
else
forj=i:
-1:
1;
tmpavg=tmpavg+rudedata（j）;
end;
tmpavg=tmpavg/i;
end
rudedata（i）=tmpavg;
end;
signal=rudedata;
%最小二乘法去线性
%procedureFIT（varDATA:
arrayofDouble;Ns:
integer）;
function[DATA]=FIT（sig）
SX=0;
SY=0;
Ns=length（sig）;
fori=1:
Ns;
SX=SX+i;
SY=SY+sig（i）;
end;
SXOSS=SX/Ns;
ST2=0;
B=0;
fori=1:
Ns;
T=i-SXOSS;
ST2=ST2+T*T;
B=B+T*sig（i）;
end;
B=B/ST2;
A=（SY-SX*B）/Ns;
fori=1:
Ns
DATA（i）=sig（i）-A-B*i;
end;
%去除背噪函数
function[signal]=FindSignal（rudedata）
Ns=length（rudedata）;
maxv=max（rudedata）;
k=0;
fori=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;
fori=2:
k;
minv=peak（i）;
forj=peakindex（i-1）+1:
peakindex（i）-1;
if（rudedata（j）<=minv）
minv=rudedata（j）;
valley（i-1）=minv;
valleyindex（i-1）=j;
end;
end;
end;
fori=2:
k-1;
maxv=max（valley（i-1）,valley（i））;
forj=valleyindex（i-1）:
valleyindex（i）;
rudedata（j）=rudedata（j）-maxv;
ifrudedata（j）<0
rudedata（j）=0;
end;
end;
end;
signal=rudedata;
%平滑谷地函数
function[signal]=Smoothness（rudedata）
Ns=length（rudedata）;
sumv=0;
k=0;
fori=2:
Ns-1;
if（rudedata（i）>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;
fori=1:
Ns;
if（abs（rudedata（i））<=avg）;
signal（i）=0;
else
signal（i）=rudedata（i）;
end;
end;
end;
附件2、连续Haar小波变换源码
functionwcoefs=myHaarcwt（
升级会员 