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

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

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

問題在於，除了打電話這種事情，人們還是有傳送數據的需求的，比如說像尋呼機這種無線數字通信系統，更早的，比如鄭君裡老師在 《教與寫的記憶-信號與系統評註》提到的他年輕時候的神器「1200波特數傳機」這東西用現在的話講叫做「1200波特率數據機」，送你一臺上網用，你肯定嫌慢，但是在當時是要國家立項的重大課題。如果你願意去一些通信原理或是信號與系統的課本裡面考古，也許會看見有些習題專門探討如何設計一個模擬的升餘弦滾降濾波器。在那個時代裡，數字電路的任務是把要發送的比特信息變換成脈衝信號，就是一些列各種幅度（多進位調製）的方波，我們在信號與系統的課程裡知道，方波信號的帶寬是無窮大的，所以後級的模擬成形濾波器負責把這些方波的頻譜帶寬變小，同時又要滿足時域採樣點無失真的準則。

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

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

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

所以，在數位訊號處理的領域，基帶成形濾波器通常的選擇還是FIR，在用數字的方法進行基帶成形時，需要額外注意的一個事情就是多速率，這是容易混淆的地方，這種成形濾波器有兩個速率，一個是輸入的符號的速率，濾波器的輸入是一系列或正或負的脈衝串，這個脈衝串是由要傳輸的信息比特經過符號映射得到的，如果是多電平調製的方式，脈衝串會有多種幅度，最簡單的映射方式是：信息比特0映射為脈衝符號-1，信息比特1映射為脈衝符號+1，這個脈衝符號是什麼呢，就是「奈奎斯特準則」裡面所說的「採樣點無失真」的採樣點，這個採樣點指的就是脈衝符號。另一個速率呢，就是濾波器輸出波形的採樣率，根據採樣定理，不能小於2倍的信號帶寬，而且，有時候為了其他的原因，比如簡化濾波器的設計，可能會用更高的採樣率。

好了，這裡我們又要展示一個蹺蹺板了，如同短時傅立葉變換中的時域、頻域解析度是蹺蹺板的兩端。這個基帶成形設計也有個蹺蹺板，讓我們來看看。

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

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

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

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

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

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

滾降係數0.5

插圖，時域波形

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

插圖，濾波器頻響

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

插圖，加窗DFT頻譜

然後我們關注一下短時的幅頻特性，上STDFT（短時傅立葉變換）頻譜。

插圖，瀑布圖頻譜

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

磚牆是不可能完美實現滴，無限逼近吧，你看這個頻譜，還記得信號與系統課本裡面，有個詞彙叫做「吉布斯現象」麼？

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

插圖，濾波器頻響

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

插圖，加窗DFT頻譜

短時分析也是如此

插圖，瀑布圖頻譜

接下來，再把滾降係數為改為1。

插圖，時域波形

插圖，濾波器頻響

插圖，加窗DFT頻譜

插圖，瀑布圖頻譜

通信系統的確是複雜且零碎的一門學科，希望本文對您的學習能盡到綿薄之力，網絡上曾有詩曰：

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

以下是代碼片段：

%///////////////////////////////////////////////////////////% 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]);