深度剖析MATLAB信號處理仿真-基帶脈沖成形

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

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

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

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

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

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

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

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

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

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

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

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

下面我們開始玩弄這個代碼了。

先找一個比較中庸的配置，滾降系數(shù)0.5，采樣率是3倍的符號率，濾波器延遲是8（個脈沖符號），由于FIR是中心對稱的，你可以看到，總共的抽頭系數(shù)是3*8*2+1=49個，通過下面這張時域的圖你可以看到濾波器沖擊響應(yīng)的時域波形，以及輸入的信息比特被映射為幅度是+1、-1的脈沖串，并且，被做了補零插值當做濾波器的輸入。然后，在濾波器輸出信號中你看看，是否實現(xiàn)了采樣點無失真呢（請關(guān)注那些幅度是+1，-1的樣點）？這個時域的信號的原始信息比特是0x47，數(shù)字衛(wèi)星標準中的著名幀頭。噢，等等，你是否覺得有點眼熟，對了，其實數(shù)字化的成形濾波器就是一個多速率信號處理里面的插值濾波器，這個我們在之前的系列文章中是有討論的。

滾降系數(shù)0.5

插圖，時域波形

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

插圖，濾波器頻響

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

插圖，加窗DFT頻譜

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

插圖，瀑布圖頻譜

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

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

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

插圖，濾波器頻響

這信號帶寬是600Hz，帶外雜散很大。

插圖，加窗DFT頻譜

短時分析也是如此

插圖，瀑布圖頻譜

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

插圖，時域波形

插圖，濾波器頻響

插圖，加窗DFT頻譜

插圖，瀑布圖頻譜

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

天若有情天亦老，人學(xué)通信死得早。商女不知亡國恨，隔江猶看信息論。問君能有幾多愁，誤碼率都不會求。兩岸猿聲啼不住，互相討論譜密度。忽如一夜春風(fēng)來，信號流圖不會排。洛陽親友如相問，直接去問韓聲棟。風(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]);