学习相关算法.docx
《学习相关算法.docx》由会员分享，可在线阅读，更多相关《学习相关算法.docx（11页珍藏版）》请在冰豆网上搜索。
学习相关算法
写个故事捧场。
上世纪90年代，BP机还很流行，电信是绝对领导地位，用户很多，同频点发射机也很多。
当时有个问题，同一频点的CALL台发射机多了以后，由于每个发射机到交换机的传输延迟不一样，导致覆盖重叠区域，到达用户BP机的信号相位差过大，误码很高。
如果是放在现在，或者当时的欧美（摩托的工程师如是说），这个问题解决得很简单也很完美，就是每个发射机那放个GPS时钟，估测个最长延迟，每个发射机加个缓存，全部在同一时刻发同一条寻呼信息就结了。
可那个年代可还是第一次海湾战争结束不久的年代，GPS是跟战斧联系在一起的，是太贵还是怎么，具体细节不记得了，总之，这个方案当时不现实。
再强调一次，年代久远，不是直接负责人，细节真不记得了，经不起推敲的不要质疑。
但原则上肯定是没错的，这个项目给我印象太深了。
这个问题讨论后方案的核心就落在如何测量每个发射机链路的延迟上，这个延迟要测得比较准，大概是不大于寻呼数据帧每一位的1/4位宽，也就是为了保证相位差了。
具体怎么测呢？
大概就是夜深人静之际，关掉所有寻呼发射机，然后逐台发一特殊数据帧，该帧从传呼中心交换机-有线传输链路-发射机-空中-接收机-鉴别电路，算出收发时间差。
电信很配合，与电信相关的几个传输环节的细节都搞得很清楚了，不会有问题；空中这一块直观上，BP都可以收到，只是重叠区域有误码，那只开一台发射机接收肯定是没问题了。
于是第一方案看上去很朴素，做一个特定码流发出去，最后无线解码出来，允许有一部分被干扰，总有没被干扰的。
事实上一个帧重复发几次就行了，是的，默认的前提就是不能全部被干扰。
找到后测脉冲个数，周期，与自己发出去的相同则认为真正找到，自然可测出整个回路的传输延迟。
第一方案看上去也没有什么大问题。
方案定了，这个东西做起来相对简单，就生成测试码流，接到传呼中心发出去，无线接收机是现成的，然后解码出来的测试码流识别电路搭下。
当然还有花里胡哨的上位机，测好所有延迟后在传呼中心建的同步缓冲器。
这个项目做容易，现场测试是个麻烦事。
第一，只能后半夜，当然电信几乎所有的割接都是后半夜。
第二要关几十台发射机（真有这么多，我也不知道为什么，印象里好多覆盖全城的小寻呼台就1，2台发射机），动静很大（细节没问），每次去都是10几个人。
第三，每晚只能中断正常业务5~10分钟。
总之，边测边调是不现实的。
第一次出师不利，这符合大家的心理预期。
有点不妙的是测试数据，跟预期的差得比较大。
由于是第一次测，局方很重视，给的时间比较多，那天晚上前后测了半个小时，结果发现几百次测试数据，只有几次能找到发出的信号，其他要不变形严重，要不干脆就淹没不见了。
恩，看来开始还是太乐观，是按BP机都可以正确解码考虑的。
我们用的是专用的寻呼检测接收机，设备应该不是问题，先想到的是改进原来的电路，但方向还是提高信噪比。
再试，检出率高了，但误检率也高了，同个发射机上10次检出数据中居然都没有一致的多数。
回来改，再去试，反复多次后，人家脸色都不好看了，这才觉得可能此路不通了。
期间请大学教授，相关厂家，都没有好的解决办法。
写到这里，我觉得这确实是个问题，平时BP机都可以正常工作的，怎么专用测试接收机反而干扰这么严重？
由于我不是那个项目组的，只是该项目组手段用尽，严重受阻，提请全体开发人员建议，最后归结到：
接收机出来就这么个波形，能不能准确找到我们发出的信号？
这才是我可以考虑的问题。
相关检测的想法并不是我提出来的，而正是他们项目组的一个原来在摩托的售后工程师。
他来我们这小公司，属于英雄难过美人关的故事，先按下不表，呵呵。
碰到这个问题后，他想起他曾经遇到过的相似的问题，用的是相关检测的方法。
他大学毕业，呆的是西安水下兵器研究所，就是研制鱼雷的地方。
80年代中期，中美蜜月期，美国佬把他第二好的鱼雷MK48所有的技术，包括细节都转让给了中国，具体接手，就是他们所。
据他说，图纸有好几吨，工艺绝佳，所有圆形成品PCB，拿铅笔顶在圆心可以保持平衡。
而且就像ZLG前几天所说的那样，所有技术从最基本的原理讲起，一直到具体实现，都很细致，这点看，美国人很实在。
。
如果是安心搞技术的人，这地方真不错，埋下头啃就是了。
可惜他不是，呆了2年去了摩托。
不过耳濡目染，他记得那鱼雷制导的方式之一就是发出超声波，接收回波后用相关检测判断距离的。
海里面有很多干扰，回波也是乱七八糟的，但是那相关器准得很哦。
去了摩托，居然还是碰到了相关检测，这就是大名鼎鼎的CDMA。
这位老兄记住了培训中那相关函数的经典图像：
从基底噪声中脱颖而出的有效信号。
却依然不知道怎么让它脱颖而出。
。
如今要用了，于是现学。
打电话给西安的师傅：
怎么做相关检测啊？
师傅：
相关检测啊，先相乘后积分啊！
同事：
先相乘后积分。
。
噢。
。
和各位看官一样，我和我那同事也是一头雾水，完全不知所云。
我不是科班，但公司里十几号搞开发的也说不清。
看这个帖子的回应也可以知道，即便是科班，不专门做这行，也会选择性遗忘的，呵呵。
接着说项目，很奇怪，已经是山穷水尽了（我没着力描写，实际拖了半年，很令人沮丧的,特别是他们项目组几个人），有人提个很有意思的思路，居然没人去落实，包括他自己，打个电话没明白似乎就算了。
我却对这个概念很感兴趣，什么方法可以检测到完全淹没到噪声里的信号？
回头翻书，发现信号与系统或数字信号处理的头1,2章就提到相关函数的概念和性质，但也几乎都是一笔带过。
不被项目弄到焦头烂额，很难让人去关注它。
我大一的高数还是蛮不错的，仅仅几年没用，看到那些+-无穷积分就晕乎乎的了，看到相关函数的表达式和一些推导完全没有感觉，无法领会那些符号背后的物理意义。
直到几天后无意中看到一本大专教材的关于CDMA原理的简单介绍。
这本书里的CDMA正交码和相关器的介绍很简单，举了最简单的一组4个正交码：
（1111,1-11-1,11-1-1,1-1-11），然后说，所谓正交码就是任意2串码的相关系数都等于0，所谓相关系数就是对应位相乘再相加：
1*1+1*（-1）+1*1+1*（-1）=0。
耳边突然响起”相关检测啊，先相乘后积分啊！
“于是突然懂了所谓正负无穷积分的相关函数表达式。
一周后，电路改进完毕，测试一次性通过。
码了一晚上，我看看能不能用简单的例子把原理说明白：
还是上面4个正交码（1111,1-11-1,11-1-1,1-1-11），其他码也一样。
其中1就是高电平，-1就是低电平，就像232电平一样。
上面说了，相关系数=对应位相乘再相加1*1+1*（-1）+1*1+1*（-1）=0，那相关函数到底是什么呢？
还是简单点，虽然不严谨，但绝对是这个意思：
假设我发了一个11-1-1，开始从接受到的码流里找，假设收到的码流跟上面排列一样。
那么，收到4个后进行第一次比较，算第一个相关系数：
f（0）=1*1+1*1+1*（-1）+1*（-1）=0；
再收一个（还是1），再比：
f
（1）=1*1+1*1+1*（-1）+1*（-1）=0；
再收一个（是-1），再比：
f
（2）=1*1+1*1+1*（-1）+（-1）*（-1）=2；
下一个
（1）：
f（3）=1*1+1*1+（-1）*（-1）+1*（-1）=2；
f（4）=1*1+（-1）*1+1*（-1）+（-1）*（-1）=0；
f（5）=（-1）*1+1*1+（-1）*（-1）+1*（-1）=0；
f（6）=1*1+（-1）*1+1*（-1）+1*（-1）=-2；
f（7）=（-1）*1+1*1+1*（-1）+（-1）*（-1）=0；
f（8）=1*1+1*1+（-1）*（-1）+（-1）*（-1）=4；
f（9）=1*1+（-1）*1+（-1）*（-1）+1*（-1）=0；
f（10）=（-1）*1+（-1）*1+1*（-1）+（-1）*（-1）=-2；
f（11）=（-1）*1+1*1+（-1）*（-1）+（-1）*（-1）=-2；
f（12）=1*1+（-1）*1+（-1）*（-1）+1*（-1）=0；
上面的步骤就算出了相关函数的值，其中算出的最大值4对应的自变量8就是特征码返回的时刻。
实际的运用，无非是特征码（11-1-1）更长一点（更合理一点，随机码最好），也不是这种1位A/D的，接收到的码流也更长一点。
当长度大到一定程度，这种算法就太笨了。
这时相关函数的计算可以化为卷积计算，卷积计算可以用FFT。
。
我前面故事中的项目就是这么傻算出来的，速度可以接受，再说，那时没人会算，效果达到了就行了。
。
85楼讲解了特征码的相关，我也稍微再讲解一下伪随机码的应用。
我不是搞语音的，也不是做测距的，做测距纯属好玩，也为了写讲义。
对于伪随机序列没有多少研究，只做一点原理介绍。
我上面提到的chirp+相关的测距方法，虽然有很多优点，但是系统的软硬件都比较复杂，不太适合低档的mcu.而一般的测距方法（如老叶前面提供的）抗干扰能力较差，提高灵敏度或有其他干扰产生一个干扰脉冲就会导致测量错误。
然而把伪随机序列应用到这种方法后，可以极大的提高抗干扰的能力。
产生一个伪随机序列，假设101100111000，用来控制发射脉冲，比如1则发送一个周波（周期为1/40K）的信号，0则停止一个周波。
回波放大解调整形后，又得到了一个序列。
把接收到的序列与伪随机序列101100111000做相关，可以排除干扰。
序列相关值的计算，teddeng已经用例子说的非常清楚，仔细的看看就会明白。
相关值c可以用下面的公式说明:
c=（a–b）/（a+b）
a是伪随机序列与比较的序列逻辑相同的数量，b是不同的数量。
因此，如果接受序列也是101100111000，即b=0,那么相关值c=1，说明完全相关;同理，如果接受序列是010*********，即a=0, 那么相关值c=-1，说明两个序列无关。
现在，我们只是把原来发送若干周波的信号然后检测上升沿，改为发射一串调制的伪随机序列信号，然后计算相关值，设定相关值大于某个值，比如0.8，则认为接受正确，这样可以排除干扰脉冲，极大的提高可靠性。
对大致LZ方法理解，使用大家容易懂的白话描述如下:
1、使用单次脉冲的方法测量受制于脉冲能量和接受灵敏度的原因，所以测量距离降低。
2、使用LZ的办法，发送一串数字序列，这个序列中1的个数将决定在这个数字序列的时间内的发送能量，在接受端接受的已经不是所谓的脉冲，而转化对这个数字序列的能量进行判断，LZ认为只要设定一个阀值，例如80%的能量系数，只要检测到返回这么多能量，就视为一次发射。
如果以上两点正确，则会存在一些显而易见的问题。
实际做法上，回波可以被转化为一个长长的序列，如...00000000001000000000....其中的1是一个干扰。
而相关值计算非常简单，而且是递推的，在mcu的定时中断里可以完成。
一但相关值大于设定，则检测峰值，峰值处可以被认为是检测到了真正有效的回波。
当然，如果对成本没有苛刻要求，用高速cpu+RAM+DA/AD发射chirp信号，然后处理数字信号是更好的方案。
//产生一个伪随机序列，假设101100111000，用来控制发射脉冲，比如1则发送一个周波（周期为
//1/40K）的信号，0则停止一个周波。
回波放大解调整形后，又得到了一个序列。
把接收到的序列与伪随机序
//列101100111000做相关，可以排除干扰。
纯理论的东西，有可行性否？
只输入一个40K频率的脉冲，对于数字滤波器，或硬件谐振滤波器都还没能反应过来吧？
粉丝是对的。
伪随机序列用来调制，一个周波，数字滤波器还能勉强做到，但对于硬件谐振滤波器估计很难，应根据解调芯片的要求确定周波数，是我简单化了。
85楼的卷积运算实例，公式为：
f（x）=b0*a0+b1*a1+b2*a2+b3*a3，
其中，a0,a1,a2,a3为发射码信号串，即1，1，-1，-1，
b0,b1,b2,b3为接收码信号串，每接收一码，信号串向前移出一位（丢弃），尾部再补入接收到的
这一位新码，再安照接收信号码（时间）顺序，依次算出f（0）,f
（1）,f
（2）,f（3）,f（4）,f（5）.......
找出这些函数值的最大值，其对应的码（时间）值，即为返回的时刻。
85楼的卷积运算实例中，上面的步骤就算出了相关函数的值，其中算出的最大值4对应的自变量8就是特征码返回的时刻。
实际上是将这一个函数看作区间的指示函数，利用卷积的“滑动平均”功能，求出最大值。
卷积和循环卷积的区别：
7个1卷积7个1，也就是1111111卷积1111111
用卷积来计算（为了对齐，这里加入了不少0）
卷积：
00000001111111
00000001111111
--------------*乘
00000001111111
0000001111111
000001111111
00001111111
0001111111
001111111
01111111
--------------+加
01234567654321
注意：
加完后不要进位,这里前一个不是后一个的十倍,而是一个平等的信号。
循环卷积：
00000001111111
00000001111111
--------------*乘
00000001111111
0000001111111
000001111111
00001111111
0001111111
001111111
01111111
-------------->>超出部分移动到右边
00000001111111
00000001111111
00000001111111
00000001111111
00000001111111
00000001111111
00000001111111
--------------+加
00000007777777
对于楼主的（a0a1a2a3）*（b0b1b2b3）
卷积：
a0 a1 a2 a3
b0 b1 b2 b3
--------------------*
a0b3a1b3a2b3a3b3
a0b2a1b2a2b2a3b2
a0b1a1b1a2b1a3b1
a0b0a1b0a2b0a3b0
-----------------------------------+
a0b0,a0b1+a1b0,a0b2+a1b1+a2b0,a0b3+a1b2+a2b1+a3b0,a1b3+a2b2+a3b1,a2b3+a3b2,a3b3
循环卷积结果：
a0b3a1b3a2b3a3b3
a1b2a2b2a3b2a0b2
a2b1a3b1a0b1a1b1
a2b0a0b0a1b0a2b0
--------------------+
a0b3+a1b2+a2b1+a2b0,a1b3+a2b2+a3b1+a0b0,a2b3+a3b2+a0b1+a1b0,a3b3+a0b2+a1b1+a2b0
显而意见，卷积和循环卷积的个数不一样，傅立叶变换和快速傅立叶变化的结果不一样。
傅立叶变换和快速傅立叶变换，对于周期信号的测量，结果是基本相近的。
因为楼主给的序列太短了，所以用的卷积，而没有用FFT。
N=8192;
M=2*N;
Fs=44100;
Ts=1/Fs;
t=0:
Ts:
（N-1）*Ts;
window=blackman（N）';
Transmit=chirp（t,500,N*Ts,5000）.*window;
Template=Transmit（length（Transmit）:
-1:
1）;
HTemplate=fft（Template,M）;
HTemplate=HTemplate（1:
N）;
rec=audiorecorder（Fs,16,1）;
record（rec,0.8）;
pause（0.2）;
wavplay（Transmit,Fs）;
pause
（1）;
stop（rec）;
Receive=getaudiodata（rec,'single'）';
Receive=Receive（floor（0.2*Fs）:
length（Receive）-1）;
sizeOfData=length（Receive）;
%Result=zeros（1,sizeOfData）;
%window=blackman（M）';
%segment=N/4;
%fori=0:
segment:
（floor（（sizeOfData-M）/segment）*segment）
% FRec=fft（Receive（（i+1）:
（i+M））.*window）;
% Corr =HTemplate.*FRec（1:
N）;
% CorrFull=[Corr,conj（Corr（N:
-1:
1））];
% Result（（i+1）:
（i+M））=Result（（i+1）:
（i+M））+（real（ifft（CorrFull）））;
%end;
Result=conv（Receive,Template）;
Result=abs（Result）;
[max1,t1]=max（Result）;
WindowSize=50;
[max2,t2]=max（Result（（t1+WindowSize）:
length（Result）））;
figure
（1）;
subplot（2,1,1）,plot（Receive）;
subplot（2,1,2）,plot（（Result））;
if（max2>max1*0.05）
d=340*（（t2+WindowSize）*Ts）/2;
title（['Distance:
',num2str（d,'%6.3f'）,'m']）;
end
大家可以下载一个matlab免安装绿色版本运行此程序。
被%注释的部分是与Result=conv（Receive,Template）作用相同的FFT/IFFT卷积。
而matlab里的conv卷积命令据称也是用FFT/IFFT实现。
这个方法对于一般运动物体也有效，只是误差可能较大，因为chirp信号时间较长，如程序中8192/44100=185ms。
对于30m/s的物体，185ms会移动5.6m,失去了测距的意义。
chirp信号是扫频信号，所以多普勒效应的影响是把chirp信号左移或右移，绝大部分还是与模板信号相关。
下面我再对卷积与FIR做一些补充，这是我答应过的。
数字滤波器可以分为FIR与IIR两种,两种滤波器的核心（kernel）都是冲击响应，两者之间的差别在于是否使用了迭代（递推）。
IIR的应用公式可以通过Laplace传递函数变换到迭代形式，一个系统的冲击响应的Laplace变换就是系统的传递函数，例如以前提到的RC滤波器的例子，其IIR迭代公式为：
y（k）=（1-Ts/τ）*y（k-1）+Ts/τ*x（k）
其中τ=RC,Ts为采样周期。
卷积
前面讲过，一个系统的输出可以通过输入与系统冲击响应的卷积得到，这是FIR滤波器计算公式核心：
y（t）=∫x（t-τ）h（τ）dτ=x（t）◙h（t）
卷积的离散表达式为：
y（k）=Σx（k–i）*h（i）=x（k）*h（0）+x（k-1）*h
（1）+…+x（k-N）*h（N）
上式中的h（i）是系统的冲击响应序列乘以采样间隔Ts。
其物理意义是，一个输入x（t）可以被分割成无数个脉冲Σx（k）,每个脉冲都会产生一个冲击响应h（t），这些冲击响应叠加在一起，就形成了连续的输出y（t）,如同顺序击打无数面鼓所产生的声音。
1.
2.第k-N个脉冲的冲击响应:
x（k-N） *[h（0） h
（1） … h（N）]
3.第k-N+1个脉冲的冲击响应:
x（k-N+1）* [h（0） h
（1） … h（N）]
4. …
5.第k-1个脉冲的冲击响应:
x（k-1） * [h（0） h
（1） … h（N）]
6.第k个脉冲的冲击响应:
x（k） * [h（0） h
（1） … h（N）]
7.第k+1个脉冲的冲击响应:
x（k+1） * [h（0） h
（1） … h（N）]
8. ↑ ↑
9. y（k） y（k+1）
复制代码
因而，在第k时刻的输出为：
y（k）=x（k）*h（0）+x（k-1）*h
（1）+….+x（k-N）*h（N）
FIR滤波器
前面讲过一些FIR的基本设计知识，比如利用设计一个窗口频谱然后用逆FFT得到FIR的系数，或者由已知的公式例如sinc函数求出FIR系数。
sinc函数是一个有趣的函数，它的频率响应是一个窗口，而一个频谱是sinc函数所对应的时域函数则是一个窗口，即
h（t）=1, t h（t）=0, t>=Tc
离散形式就是：
y（k）=x（k）+x（k-1）+….+x（k-N-1）=[1,1,…,1]*X
这个公式很熟悉吧？
这就是移动平均算法。
对移动平均算法的内核序列[1,1,…,1]/N作FFT，可以求得频谱：
可以看出，移动平均算法的频谱是一个sinc函数，移动平均滤波器是一个FIR低通滤波器。
阻容FIR与IIR滤波器
RC电路对于大家并不陌生，其阶越响应，通俗地说，直流充电过程是以一个电容电压（假设充电电压为1V）逐渐接近直流电压的过程，用数学公式表示为：
升级会员 