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