深度剖析MATLAB信號(hào)處理仿真-基帶脈沖成形

本次我們探討另外一個(gè)在本科階段讓我們頭痛的東西，通信原理之必考曲目，拼死也要背下來的內(nèi)容，基帶脈沖成形。然而俺對(duì)這個(gè)東西的理解和認(rèn)識(shí)卻是在本科以后的事情。

早年（比如摩托羅拉手機(jī)時(shí)代）的基帶成形都是用模擬電路做的，那會(huì)兒的數(shù)字電路密度極低，想想大學(xué)本科數(shù)電實(shí)驗(yàn)里面的各種74系列芯片，如果用這個(gè)東西拼個(gè)數(shù)字濾波器估計(jì)會(huì)瘋掉。況且，就算有數(shù)字濾波器，高速高精度的ADC、DAC也是個(gè)問題。所以，早些年的數(shù)字電路課本的名字通常叫做“脈沖與數(shù)字電路”，言下之意，這東西用來處理脈沖信號(hào)的，而且，也就處理處理脈沖信號(hào)，千萬別想著整太復(fù)雜的東西，那會(huì)還是一個(gè)模擬電路統(tǒng)治著通信系統(tǒng)的時(shí)代。

問題在于，除了打電話這種事情，人們還是有傳送數(shù)據(jù)的需求的，比如說像尋呼機(jī)這種無線數(shù)字通信系統(tǒng)，更早的，比如鄭君里老師在 《教與寫的記憶-信號(hào)與系統(tǒng)評(píng)注》提到的他年輕時(shí)候的神器“1200波特?cái)?shù)傳機(jī)”這東西用現(xiàn)在的話講叫做“1200波特率調(diào)制解調(diào)器”，送你一臺(tái)上網(wǎng)用，你肯定嫌慢，但是在當(dāng)時(shí)是要國家立項(xiàng)的重大課題。如果你愿意去一些通信原理或是信號(hào)與系統(tǒng)的課本里面考古，也許會(huì)看見有些習(xí)題專門探討如何設(shè)計(jì)一個(gè)模擬的升余弦滾降濾波器。在那個(gè)時(shí)代里，數(shù)字電路的任務(wù)是把要發(fā)送的比特信息變換成脈沖信號(hào)，就是一些列各種幅度（多進(jìn)制調(diào)制）的方波，我們?cè)谛盘?hào)與系統(tǒng)的課程里知道，方波信號(hào)的帶寬是無窮大的，所以后級(jí)的模擬成形濾波器負(fù)責(zé)把這些方波的頻譜帶寬變小，同時(shí)又要滿足時(shí)域采樣點(diǎn)無失真的準(zhǔn)則。

我們?cè)?a target="_blank">數(shù)字信號(hào)處理課程里面學(xué)習(xí)過IIR濾波器，而且還有“雙線性變換法”，“沖擊響應(yīng)不變法"，以及各種讓我們頭暈的東西，我小時(shí)候第一次學(xué)這東西的時(shí)候在想，整這個(gè)玩意兒干嘛，后來才明白，這東西是為了用數(shù)字的方法來實(shí)現(xiàn)以前的模擬濾波器，模擬濾波器都是有極點(diǎn)的，映射到數(shù)字域中，就是IIR濾波器，那么為什么要替換掉模擬濾波器呢，有兩個(gè)原因，一是為了提高通信產(chǎn)品的一致性，模擬元件比如電容、電阻的值是無法嚴(yán)格準(zhǔn)確生產(chǎn)的，至于電感就是個(gè)更加不靠譜的東西，這就導(dǎo)致每個(gè)模擬電路元件被制造出來都不會(huì)完全一樣，如果系統(tǒng)中大量使用模擬元件，最后誤差的積累會(huì)是個(gè)大問題，所以設(shè)計(jì)者會(huì)預(yù)先在電路里面加入一些參數(shù)可調(diào)的元件，在最后的組裝階段讓工人師傅看著測(cè)試儀器手動(dòng)調(diào)節(jié)產(chǎn)品。數(shù)字電路則不同，一旦被生產(chǎn)出來，行為完全一致，后期的組裝調(diào)試環(huán)節(jié)大為簡化，所以更加有生產(chǎn)的效率，若你玩過紅色警報(bào)、星際爭霸應(yīng)當(dāng)知道，快速生產(chǎn)可是個(gè)大事情。

以上這個(gè)行為一致性的優(yōu)點(diǎn)只是一方面，另一方面，更加有吸引力的是，數(shù)字電路中有一個(gè)概念叫做”等比縮小“，比如最早的計(jì)算機(jī)有房子那么大，可能連你手機(jī)里面處理器計(jì)算能力的百分之一都沒有，我以前為了做某些算法，去70年代的IEEE trans 里面考古，發(fā)現(xiàn)某些科學(xué)家在文章的末尾貼試驗(yàn)數(shù)據(jù)的章節(jié)里炫富，”俺這個(gè)快速算法在擁有1M字節(jié)內(nèi)存的XX型號(hào)的IBM豪華大型機(jī)上的運(yùn)行時(shí)間是12秒“，俺當(dāng)時(shí)感到唏噓不已。所以，你在現(xiàn)在的港產(chǎn)警匪片里面再也看不到大佬們要專門配備一個(gè)小弟替他拿電話并且大佬還要用移動(dòng)電話砸人了。

如果不是為了替換模擬時(shí)代的濾波器，或是為了適應(yīng)一些特殊的自適應(yīng)濾波器的需求，我是會(huì)盡量回避使用IIR濾波器的，尤其是在用FPGA作為實(shí)現(xiàn)手段的場(chǎng)合。為什么？因?yàn)檫@個(gè)有極點(diǎn)的濾波器真是對(duì)量化噪聲極其敏感啊，動(dòng)不動(dòng)就自激振蕩啊，調(diào)過一次之后被折騰的七葷八素的，后來就再也不敢用了，還是乖乖的用FIR濾波器，讓領(lǐng)導(dǎo)給買個(gè)大芯片頂上先。如果是在處理器上，尤其是支持浮點(diǎn)的CPU，做做倒是無妨，但是也要注意濾波器如果階數(shù)高了，IIR這東西仍然是個(gè)危險(xiǎn)品。

所以，在數(shù)字信號(hào)處理的領(lǐng)域，基帶成形濾波器通常的選擇還是FIR，在用數(shù)字的方法進(jìn)行基帶成形時(shí)，需要額外注意的一個(gè)事情就是多速率，這是容易混淆的地方，這種成形濾波器有兩個(gè)速率，一個(gè)是輸入的符號(hào)的速率，濾波器的輸入是一系列或正或負(fù)的脈沖串，這個(gè)脈沖串是由要傳輸?shù)男畔⒈忍亟?jīng)過符號(hào)映射得到的，如果是多電平調(diào)制的方式，脈沖串會(huì)有多種幅度，最簡單的映射方式是：信息比特0映射為脈沖符號(hào)-1，信息比特1映射為脈沖符號(hào)+1，這個(gè)脈沖符號(hào)是什么呢，就是“奈奎斯特準(zhǔn)則”里面所說的“采樣點(diǎn)無失真”的采樣點(diǎn)，這個(gè)采樣點(diǎn)指的就是脈沖符號(hào)。另一個(gè)速率呢，就是濾波器輸出波形的采樣率，根據(jù)采樣定理，不能小于2倍的信號(hào)帶寬，而且，有時(shí)候?yàn)榱似渌脑?，比如簡化濾波器的設(shè)計(jì)，可能會(huì)用更高的采樣率。

好了，這里我們又要展示一個(gè)蹺蹺板了，如同短時(shí)傅里葉變換中的時(shí)域、頻域分辨率是蹺蹺板的兩端。這個(gè)基帶成形設(shè)計(jì)也有個(gè)蹺蹺板，讓我們來看看。

所謂蹺蹺板，就是說，有若干個(gè)參數(shù)，而且相互之間有制約，調(diào)節(jié)了一個(gè)參數(shù)，另外的參數(shù)會(huì)隨之改變，通常我們的做法是，先固定一部分，然后調(diào)節(jié)另外的，再觀察調(diào)節(jié)的結(jié)果。

幸好，這里只有三個(gè)參數(shù)，根據(jù)通信原理的課本，當(dāng)采用升余弦滾降濾波器進(jìn)行成形的時(shí)候，這三個(gè)參數(shù)為：脈沖符號(hào)的速率 Rs，滾降系數(shù)beta，生成的窄帶信號(hào)的帶寬BW，課本上說了，BW=Rs(1+beta)/2，beta是一個(gè)介于0、1之間的小數(shù)，這個(gè)表達(dá)式很簡單，它說明了這樣一個(gè)事實(shí)：對(duì)于符號(hào)速率為Rs的脈沖串，如果要做到采樣點(diǎn)無失真，用一個(gè)帶限信號(hào)傳輸它，那么這個(gè)帶限信號(hào)的帶寬將會(huì)介于0.5Rs和Rs之間。這時(shí)我們蹺蹺板的一頭是信號(hào)帶寬，另一頭是濾波器的頻響形狀，你看下面這個(gè)圖，從wiki拷貝下來的，這個(gè)beta=0的時(shí)候是“磚墻”牌的理想低通啊，這可不是鬧著玩的，所以說，如果要求傳輸固定速率的脈沖符號(hào)，要想節(jié)省帶寬，就要做一個(gè)能夠逼近“磚墻”的高性能濾波器，否則就得多占用帶寬，這個(gè)被固定下來的符號(hào)速率呢，就是蹺蹺板的支點(diǎn)。另外，我們也可以舉出另外的案例，比如說，固定下來傳輸帶寬，這時(shí)蹺蹺板的兩頭分別是濾波器滾降系數(shù)和符號(hào)速率，這和之前的案例本質(zhì)上都是一樣的，不失一般性，我們下面把符號(hào)速率當(dāng)做蹺蹺板支點(diǎn)來討論。

此處，我們提供了一段代碼，這段代碼表達(dá)了一個(gè)1200波特/秒的脈沖串，在通過不同配置的成形濾波器后所產(chǎn)生信號(hào)的時(shí)頻特性。

成形濾波器采用的matlab函數(shù)是rcosine，除了要配置符號(hào)率、滾降系數(shù)、輸出波形的采樣率之外，還需要配置一個(gè)濾波器延遲，這個(gè)參數(shù)是干啥的？通信原理的課本沒說，答案在數(shù)字信號(hào)處理的課本中，因?yàn)槲覀冎?，F(xiàn)IR濾波器在概念上是個(gè)多抽頭的移位寄存器，問題在于一旦你通過滾降系數(shù)指定了濾波器的頻響曲線形狀，另外還要決定的是，你準(zhǔn)備用多少個(gè)抽頭的FIR來實(shí)現(xiàn)這個(gè)濾波器，抽頭越多，實(shí)現(xiàn)的FIR的頻響阻帶抑制越好，但是代價(jià)是，運(yùn)算量和延遲。OK，是的，你又發(fā)現(xiàn)有蹺蹺板了。注意，這個(gè)函數(shù)的延遲參數(shù)的單位是脈沖符號(hào)的個(gè)數(shù)，即相對(duì)于輸入端的延遲。

下面我們開始玩弄這個(gè)代碼了。

先找一個(gè)比較中庸的配置，滾降系數(shù)0.5，采樣率是3倍的符號(hào)率，濾波器延遲是8（個(gè)脈沖符號(hào)），由于FIR是中心對(duì)稱的，你可以看到，總共的抽頭系數(shù)是3*8*2+1=49個(gè)，通過下面這張時(shí)域的圖你可以看到濾波器沖擊響應(yīng)的時(shí)域波形，以及輸入的信息比特被映射為幅度是+1、-1的脈沖串，并且，被做了補(bǔ)零插值當(dāng)做濾波器的輸入。然后，在濾波器輸出信號(hào)中你看看，是否實(shí)現(xiàn)了采樣點(diǎn)無失真呢（請(qǐng)關(guān)注那些幅度是+1，-1的樣點(diǎn)）？這個(gè)時(shí)域的信號(hào)的原始信息比特是0x47，數(shù)字衛(wèi)星標(biāo)準(zhǔn)中的著名幀頭。噢，等等，你是否覺得有點(diǎn)眼熟，對(duì)了，其實(shí)數(shù)字化的成形濾波器就是一個(gè)多速率信號(hào)處理里面的插值濾波器，這個(gè)我們?cè)谥暗南盗形恼轮惺怯杏懻摰摹?/p>

滾降系數(shù)0.5

插圖，時(shí)域波形

然后我們看看頻域的幅頻曲線，分別是線性尺度和dB尺度，我們看到這個(gè)濾波器的設(shè)計(jì)截止頻率是900Hz，這符合公式的計(jì)算，而實(shí)現(xiàn)出來的濾波器的阻帶衰減是-60dB左右。

插圖，濾波器頻響

接下來，看看全部信號(hào)的加窗DFT頻譜，同樣，可以看到，在900Hz的位置衰減到-60dB

插圖，加窗DFT頻譜

然后我們關(guān)注一下短時(shí)的幅頻特性，上STDFT（短時(shí)傅里葉變換）頻譜。

插圖，瀑布圖頻譜

好，其它參數(shù)不變，我們改改滾降系數(shù)，首先是滾降系數(shù)為0的插圖，時(shí)域波形

磚墻是不可能完美實(shí)現(xiàn)滴，無限逼近吧，你看這個(gè)頻譜，還記得信號(hào)與系統(tǒng)課本里面，有個(gè)詞匯叫做“吉布斯現(xiàn)象”么？

注意看看濾波器的通帶寬度和阻帶衰減，慘不忍睹的-20dB啊

插圖，濾波器頻響

這信號(hào)帶寬是600Hz，帶外雜散很大。

插圖，加窗DFT頻譜

短時(shí)分析也是如此

插圖，瀑布圖頻譜

接下來，再把滾降系數(shù)為改為1。

插圖，時(shí)域波形

插圖，濾波器頻響

插圖，加窗DFT頻譜

插圖，瀑布圖頻譜

通信系統(tǒng)的確是復(fù)雜且零碎的一門學(xué)科，希望本文對(duì)您的學(xué)習(xí)能盡到綿薄之力，網(wǎng)絡(luò)上曾有詩曰：

天若有情天亦老，人學(xué)通信死得早。商女不知亡國恨，隔江猶看信息論。問君能有幾多愁，誤碼率都不會(huì)求。兩岸猿聲啼不住，互相討論譜密度。忽如一夜春風(fēng)來，信號(hào)流圖不會(huì)排。洛陽親友如相問，直接去問韓聲棟。風(fēng)蕭蕭兮易水寒，調(diào)制解調(diào)各種難。垂死病中驚坐起，明天要考載噪比！

以下是代碼片段：

%///////////////////////////////////////////////////////////% base band shaping use raised cosine FIR filter%///////////////////////////////////////////////////////////close all;clear;clc;

rate_symbol = 1200 ; % input symbol raten_itp = 3 ; % num of interpolationfilter_delay = 8 ;roll_off = 1.0 ;

head_bit = [0 1 0 0 0 1 1 1 ].';total_bit_len = 2048;bit_map = [-1 1] ;

head_bit_len = length(head_bit);info_bit = randint(total_bit_len-head_bit_len,1);

total_bit = [head_bit; info_bit];

head_itp_len = head_bit_len * n_itp;

fs_filter_out = n_itp*rate_symbol % filter output sample rate

filter_delay_in_sample = n_itp*filter_delay ;

coeff = rcosine(rate_symbol, fs_filter_out, ... 'fir/normal', roll_off, filter_delay);

coeff = coeff .' ;

% map info_bit to multi-level signalsignal_1x = bit_map(total_bit + 1) ;

len_1x = length(signal_1x) ;

len_itp = len_1x * n_itp ;signal_itp = zeros(len_itp,1) ;% write the original signal value into itp signalsignal_itp(1:n_itp:len_itp-n_itp+1) = signal_1x;

data_conv = conv(coeff, signal_itp);data_conv = filter(coeff,1, signal_itp);filter_out = data_conv ;

signal_itp_x_axis = [1:head_itp_len] + filter_delay*n_itp;

figure;head_itp = signal_itp(1:head_itp_len);head_conv = data_conv(1:head_itp_len+filter_delay_in_sample);

subplot(3,1,1); plot(coeff, 'linewidth',1.5); grid on; ylim([-0.3, 1.2]);title('Filter Coeff', 'FontSize', 16);

subplot(3,1,2); stem(head_itp,...'--ms', 'MarkerEdgeColor','k','MarkerFaceColor','g', 'MarkerSize',6 );title('Header Bits Mapped to Symbol with Interpolation', 'FontSize', 16);ylim([-1.2, 1.2]);grid on;

subplot(3,1,3); stem(signal_itp_x_axis, signal_itp(1:head_itp_len), ...'--ms', 'MarkerEdgeColor','k','MarkerFaceColor','g', 'MarkerSize',6 );grid on; hold;subplot(3,1,3); plot(head_conv, 'linewidth',1.5);grid on;title('Header Bits after Pulse Shaping Filter Output ', 'FontSize', 16);

kaiser_beta = 8 ;

kaiser_win_spectrum_plot(fs_filter_out, filter_out, kaiser_beta);title('Total Signal Spectrum', 'FontSize', 16);

len_w = floor(total_bit_len /8);len_o = floor(len_w * 0.8);len_d = len_w ;figure ;% use tic toc get running timespectrogram(filter_out, len_w, len_o, len_d, fs_filter_out);title('Spectrogram', 'FontSize', 16);% get filter frequency reponse vector over 0~pi[f_rp_vec ,w_pi] = freqz(coeff);x_half_fs = w_pi/pi(fs_filter_out/2);f_rp_vec_norm_dB = 20log10(abs(f_rp_vec)+1E-10);f_rp_vec_norm_dB = f_rp_vec_norm_dB - max(f_rp_vec_norm_dB);

figure;ideal_f_rp_vec = (x_half_fs < fs_filter_out/(2n_itp))n_itp;

subplot(2,1,1);h1 = plot(x_half_fs, abs(f_rp_vec), ...x_half_fs, ideal_f_rp_vec, '-r', 'linewidth', 1.5);grid on;

title('Filter frequency response, Linear Scale', 'FontSize', 16);legend('Actual filter response', 'Ideal filter response');subplot(2,1,2);plot(x_half_fs, f_rp_vec_norm_dB, 'linewidth', 1.5);grid on;title('Filter frequency response, dB Scale', 'FontSize', 16);ylim([-120, 5]);