Try   HackMD

通訊實驗matched filter結報

tags: 通訊 matlab
撰寫時間 : 2022/05/13

實驗目的

  • 接收端利用匹配濾波器,配合適當取樣時間來解調傳送端傳送的的兩波形所代表的二位元信號。
  • 如何模擬實現匹配濾波器。

matched filter 實驗原理

  1. 首先binary sequence
    bk{0,1}
    ,會進行訊號調變(modulation),這裡使用pulse amplitude modulation(PAM),首先每個bit通過symbol mapper,將每個bit映射到一個數值,分為以下兩種
    In antipodal signaling,ak{1,bk=11,bk=0ak=2bk1In on-off signaling,ak{1,bk=10,bk=0ak=bk
  2. 接下來通過一個pulse shaping filter
    g(t)
    ,一般
    g(t)
    使用方波,調變後的訊號
    s(t)
    就是
    akg(tkT)
    ak
    就是代表這個波型的振幅。
  3. 通過一個channel,一般假設是AWGN channel,因此接收端接收到的訊號就是
    r(t)=s(t)+w(t)
  4. 接收端通過一個filter,為了得到極大化SNR值,這個filter會根據輸入訊號
    s(t)
    的波型,設計一個match
    s(t)
    的波型,設計準則為將輸入訊號先做反轉再做延遲
    g(t)=s(inverse(tTdelay T))
  5. 取樣一段時間後才會輸出訊號,當取樣時間
    t=T
    ,此時matched filter等於correlator,而若
    g(t)
    又是方波,matched filter等效上積分器,而從電子學可知積分器可以過濾高頻雜訊。
  6. 最後是一個判斷機制,判斷這個數值是0和1,需要設計一個threshold
    γ
    ,極小化位元錯誤率
    Pe
  • 最大化
    Y(T)
    (當
    Y(t)
    取樣時間
    t=T
    )的SNR,使對抗雜訊的能力增加

    首先先正規化received filter的限制為
    0T|g(t)|2dt=T
    Y=g(t)r(t)=0Tg(Tτ)[s(τ)+w(τ)]dτ=0Tg(Tτ)s(τ)dτdesird part(signal)+0Tg(Tτ)w(τ)dτnoise partgoal :max{SNR}=max{PsignalPnoise}
    要求出noise的power
    Pnoise
    就是取其變異數(二階動差 - 一階動差平方),代表是AC的power(總功率的power - DC的power)。
    Let N=0Tg(Tτ)w(τ)dτmean of N:E[N]=0Tg(Tτ)E[w(τ)]=0dτ=0variance of N:Var{N}=E[N2](E[N]=0)2=RN(0)=F1{SN(f)}=SN(f)e+j2πf0df=SN(f)df=|G(f)|2SW(f)dfWSS process pass LTI system G(f)=N02|G(f)|2df=N02|g(t)|2dtParseval's theorem=N020T|g(t)|2dtconstraint term0T|g(t)|2dt=T=N020T|g(Tt)|2dtchange of variables
    因此SNR可以改寫為
    SNRY=PsignalPnoise=[0Tg(Tτ)s(τ)dτ]2N020T|g(Tt)|2dt
    由於底下該項是常數,因此使用柯西不等式,可以求出SNR所能到達最大值為
    SNRY=[0Tg(Tτ)s(τ)dτ]2N020T|g(Tt)|2dt[0T|g(Tτ)|2dτ][0T|s(τ)|2dτ]N020T|g(Tt)|2dt=2N00T|s(τ)|2dτonly depends on energy of s(t)
    等號成立的條件為兩向量平行,即為
    g(Tτ)=cs(τ)
    或是
    g(τ)=cs(Tτ)=cs((τT))
    ,物理意義就是將輸入訊號先做反轉再做延遲,此received filter設計就是matched filter。

lab1

問題 : duration設定為4 symbol periods, 畫出在

EbN0=10dB
EbN0=3dB
的transmitted waveform和received waveform並比較差異。

EbN0=10dB
EbN0=3dB
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

定義信噪比

SNRoPsignalPnoiseEbN0 ,其中
N0
代表高斯雜訊的一個定值(
Sw(f)N02
)、
Eb
代表on-off keying每個bit的平均能量
Eb=12A2T+120=A2T2
,因此降低SNR,
N0
保持一致,等同降低降低每個bit傳送時的平均能量,由於
Pe=Q(SNR)
,因此BER會變差。

而範例程式這邊與我上述直覺理解不同點是在維持

Eb下,調動noise的power
σ2
,降低SNR,會讓
σ2
變大,因此接收端(RX)訊號右圖會比左圖變化幅度還大。

lab2-1

問題 : 解釋程式Data_bit = (rand(1, data_number) > 0.5 )為何可產生random bits?bit 0 和bit 1,其發生機率理論上應為何?實際於模擬是否符合理論預測?

指令rand(M,N)代表返回一個

M×N的矩陣,每個element的值是取區間0到1的均勻分布$ U \sim (0, 1)$並使用logic判斷式判斷是否大於0.5,若符合回傳1,反之回傳0,回傳的data type是logic

由於是均勻分布,0到0.5的區間機率是

12,判定為logic 0;0.5到1的區間機率是
12
,判定為logic 1,實際模擬如下圖
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

藉由做很多次的隨機試驗,可以得到relative frequency
fn(A)=Nn(A)n
,當試驗次數越多,會趨於機率,記為
limnfn(A)=P(A)
,由上圖可知data number越多,機率會趨於
12

lab2-2

問題 :

data_number = 10^5; % # of bits
Fs = 10; % sampling frequency (used to generate received samples)
Data_bit = (rand(1,data_number) > 0.5 ); % random bits
p1 = ones(1,Fs); % discrete‐time rectangular pulse that represents one symbol
Data_pulse_array = (p1') * Data_bit;
Data_pulse = reshape(Data_pulse_array,1,length(p1)*data_number);

其中Data_pulse = reshape(Data_pulse_array, 1, length(p1)*data_number)代表意義。

Data_bit為由機率各半(

12)的logic 0和1形成
1×105
的矩陣,Data_pulse_array為將Data_bit的值拷貝成
10×105
的矩陣,此時data type是double,最後reshape
10×105
的矩陣,以column為順序,轉換為(1,length(p1)*data_number)的矩陣,代表意義是每10個單位傳送一個symbol。

lab3

問題 : 沒有雜訊下,matched filter output

y(T)理論值與模擬值為何?

Data_receive = Data_pulse; % received samples
D_filtered = conv(Data_receive,p1); % MF output
D_demapping = D_filtered(10:10:end) / 10; % sampling at symbol rate

首先計算兩方波做convolution,如下圖所示

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

兩個方波做convolution等於一個三角波,數學式記為
[AΠ(t2W)][AΠ(t2W)]=A22WΛ(t2W)
由於本題假設沒有雜訊,因此傳送端訊號等於接收端訊號,接下來與matched filter做convolution計算為根據上述知識,使用
A=1,2W=10
,可計算結果,使用matlab模擬驗證,假設傳送bit為1101,matched filter output
y(t)
如下
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

y(t)
t=T
時刻的時間,即上圖中simple index在10、20、30、40下的值,由於大小是
A22W
,因此需要除以一個scaling
A2W
轉為原本波型的大小
A
,本題就是
A2W=125=10
,最後繪製如下圖,藍點代表
y(T)
是時間在10、20、30、40下的數值除以一個10的scaling。
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

lab4

問題 : 解釋D_demapping=D_filtered(10:10:end)/10

同lab3問題,取matched filter在

t=T下的輸出,並除上一個
A2W=125=10
的scaling。

lab5

問題 : 解釋sgma = sqrt( 0.5/EbN0/2*10 )

如lab1觀念,定義信噪比

SNRoPsignalPnoiseEbN0 ,其中
N0
代表高斯雜訊的一個定值(
Sw(f)N02
)、
Eb
代表on-off keying每個bit的平均能量
Eb=12A2T+120=A2T2
,由於本題條件
A=1,T=1
,因此
Eb=12
,另外假設取樣頻率
Ts=10
,每個sample point的noise power為
σs2
,計算
σs2
值如下
Let N=0Tg(Tτ)w(τ)dτmean of N:E[N]=0Tg(Tτ)E[w(τ)]=0dτ=0σs2=1Tsσ2=1Ts(E[NN]E[N]2)=1Ts0T0TE[w(t)w(τ)]dtdτE[N]=0=1Ts0T0TRw(tτ)dtdτnoise is SSSWSS=1Ts0T0TN02δ(tτ)dtdτ=1TsN02N0=2Tsσs2

SNR可改寫為

SNRoPsignalPnoiseEbN0=A2T22Tsσs2=12210σs2

lab6

問題 : 畫出OOK中BER對SNR理論值與實際值做圖。

在機率equally likely

p=q=12條件下,OOK的bit error rate(BER)理論值計算為
Pe=AT2σ12πey22dy=Q(AT2σ)=Q(AT2N0T2)=Q(A2T21N0)=Q(EbN0)
與實際用matlab模擬值

 D_demap_N = (D_demapping > 0.5); % >0.5: 1; <=0.5: 0
 % BER computation
 Error_num = sum(xor(D_demap_N,Data_bit)); % same -> 0; diff -> 1
 BER = Error_num/data_number;

繪圖如下

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

由上圖可知理論值計算與matlab模擬值之間誤差非常小。

lab7

問題 : 在

EbN0=10dB
EbN0=3dB
的情況下,畫出matched filter output
y(T)
的conditional pdf
f(y(T)|bit 0 sent)
f(y(T)|bit 1 sent)
的理論值與模擬值。

index_0 = 1;
for i = 1:data_number
    if Data_bit(i) == 0 % TX == 0
        D_demapping_0(index_0) = D_demapping(i);
        index_0 = index_0 + 1;
    end
end

如果一開始傳送訊號是0,提取輸出訊號

y(T)並存入另一個自定義的array,代表選取指定事前機率下的輸出數據

histogram(D_demapping_0, 'FaceColor', '#0072BD','Normalization','pdf');

下指令histogram(X,'Normalization','pdf')畫各個數值區間內的統計數量,後面那2個參數是取數據X的probability density function。

plot(xaxis,normpdf(xaxis,0,sgma / sqrt(10)), 'LineWidth',3);

理論值計算下指令normpdf,代表使用常態(高斯)分布

N(μ,σ2),其中mean
μ
是0,而standard deviation(標準差)
σ
是需要乘以
110
的scaling,原因是根據lab3觀念,經過convolution後的訊號要除以
10
的scaling才是最終的
y(T)
y(T)
其中也包括雜訊
N(0,σ2)
,因此
σ2σ210σσ10

分別就不同SNR

EbN0下畫出conditional pdf
f(y(T)|bit 0 sent)
f(y(T)|bit 1 sent)
的理論值與模擬值。

EbN0=10dB
EbN0=3dB
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

降低SNR,在維持

Eb的條件下,調動noise的power
σ2
,會讓
σ2
變大,因此右圖的PDF會比左圖的PDF還矮胖。

程式

  • lab1
% matched filter
% ‐ on‐off keying (OOK) using rectangular pulse (T=1)
% That is, when OFF, assuming A=0, thus E1 = 0
% when ON: assuming A=1, thus E2 = 1
% The "average bit energy" (Eb) = (E1+E2)/2 = 1/2
%
clear,clc,close all;
%% parameters
data_number = 4; % # of bits
Fs = 10; % sampling frequency (used to generate received samples)

%% transmitter
Data_bit = (rand(1,data_number) > 0.5 ); % random bits
p1 = ones(1,Fs); % discrete‐time rectangular pulse that represents one symbol
Data_pulse_array = (p1')*Data_bit;
Data_pulse = reshape(Data_pulse_array,1,length(p1)*data_number);

%% AWGN channel and receiver
% AWGN channel
EbN0dB = 10; % 3 10
[a, b] = size(Data_pulse);
EbN0 = 10^(EbN0dB/10); % EbN0 is now in linear scale
sgma = sqrt( 0.5/EbN0/2*10);
noise = normrnd(0, sgma, a, b);
Data_receive = Data_pulse + noise; % received samplesa

%% generate plots
figure;
xais = 1:40;
hold on;
plot(xais, Data_pulse, 'LineWidth',2);
plot(xais, Data_receive, '--r', 'LineWidth',2);
ylim([-5 5]);
xlim([1 40]);
xlabel('sample index');
ylabel('sample value');
legend('TX waveform', 'RX waveform');
grid on;
title('waveform of OOK with Rectangular Pulse ($\frac{R_b}{N_0} = 10 \;\rm dB$)','Interpreter','latex','FontSize',15);
  • lab2
clear,clc,close all;

pro_1 = zeros(1, 20); % preallocate memory
for k = 1:10^4
    sum_1 = 0;
    data_number = k; % # of bits
    Data_bit = (rand(1,data_number) > 0.5 ); % random bits
    for i = 1:data_number
        if Data_bit(i) == 1
            sum_1 = sum_1 + 1;
        end
    end
    pro_1(k) = sum_1 ./ data_number;
end

xais = 1:10^4;
plot(xais, pro_1);
ylim([0, 1]);
xlabel('data number');
ylabel('probibilty of bit 1');
grid on;
title('probibilty of bit 1 with data number increasing','FontSize',15);
  • lab3
% matched filter
% ‐ on‐off keying (OOK) using rectangular pulse (T=1)
% That is, when OFF, assuming A=0, thus E1 = 0
% when ON: assuming A=1, thus E2 = 1
% The "average bit energy" (Eb) = (E1+E2)/2 = 1/2
%
clear,clc,close all;
%% parameters
data_number = 4; % # of bits
EbN0dB_vec = [0:2:10 11.5]; % Eb/N0 in dB
Fs = 10; % sampling frequency (used to generate received samples)

%% transmitter
Data_bit = (rand(1,data_number) > 0.5 ); % random bits
p1 = ones(1,Fs); % discrete‐time rectangular pulse that represents one symbol
Data_pulse_array = (p1')*Data_bit;
Data_pulse = reshape(Data_pulse_array,1,length(p1)*data_number);

%% no AWGN channel and receiver
Data_receive = Data_pulse; % received samples
% receiver
D_filtered = conv(Data_receive,p1); % MF output
D_demapping = D_filtered(10:10:end) / 10; % sampling at symbol rate
% decsion based on D_demapping
D_demap_N = (D_demapping > 0.5); % >0.5: 1; <=0.5: 0

%% generate plots
figure;
xaxis = 1:40;
hold on;
plot(xaxis, Data_pulse, 'LineWidth',2);
plot(xaxis, Data_receive, '--r', 'LineWidth',2);

for i = 1:4
    scatter(10*i, D_demapping(i),'b', 'LineWidth',3);
end

ylim([-5 5]);
xlim([1 40]);
xlabel('sample index');
ylabel('sample value');
legend('TX waveform', 'RX waveform', 'matched filter output at t = T','FontSize',15);
grid on;
title('matched filter output at $t = T$ with no noise channel','Interpreter','latex','FontSize',15);
hold off;

% xaxis_49 = 1:49;
% plot(xaxis_49, D_filtered, 'LineWidth',2);
% xlabel('sample index');
% ylabel('sample value');
% title('matched filter output $y(t)$','Interpreter','latex','FontSize',15);
  • lab6
% matched filter
% ‐ onoff keying (OOK) using rectangular pulse (T=1)
% That is, when OFF, assuming A=0, thus E1 = 0
% when ON: assuming A=1, thus E2 = 1
% The "average bit energy" (Eb) = (E1+E2)/2 = 1/2
%
clear,clc,close all;
%% parameters
data_number = 10^5; % # of bits
EbN0dB_vec = [0:2:10 11.5]; % Eb/N0 in dB
Fs = 10; % sampling frequency (used to generate received samples)

%% transmitter
Data_bit = (rand(1,data_number) > 0.5 ); % random bits
p1 = ones(1,Fs); % discrete‐time rectangular pulse that represents one symbol
Data_pulse_array = (p1')*Data_bit;
Data_pulse = reshape(Data_pulse_array,1,length(p1)*data_number);

%% AWGN channel and receiver
BER = zeros(1, length(EbN0dB_vec)); % pre-allocate
BER_theory = zeros(1, length(EbN0dB_vec)); % pre-allocate
for kk=1:length(EbN0dB_vec)
    % AWGN channel
    [a, b] = size(Data_pulse);
    EbN0dB = EbN0dB_vec(kk);
    EbN0 = 10^(EbN0dB/10); % EbN0 is now in linear scale
    sgma = sqrt( 0.5/EbN0/2*10);
    noise = normrnd(0, sgma, a, b);
    Data_receive = Data_pulse + noise; % received samples
    % receiver
    D_filtered = conv(Data_receive,p1); % MF output
    D_demapping = D_filtered(10:10:end) / 10; % sampling at symbol rate
    % decsion based on D_demapping
    D_demap_N = (D_demapping > 0.5); % >0.5: 1; <=0.5: 0
    % BER computation
    Error_num = sum(xor(D_demap_N,Data_bit)); % same -> 0; diff -> 1
    BER(kk) = Error_num/data_number;
    % fprintf('EbN0 in dB is %g\n',EbN0dB);
    % fprintf('Bit error rate is %g\n',BER);
    BER_theory(kk) = qfunc(sqrt(EbN0));
end

%% generate plots
figure;
semilogy(EbN0dB_vec, BER, EbN0dB_vec, BER_theory, 'LineWidth',2);
xlabel('Eb/N0 (dB), where Eb: average energy per bit');
ylabel('Bit Error Rate');
legend('OOK (Simulation)', 'OOK (Theory)','FontSize',10);
grid on;
title('OOK with Rectangular Pulse','FontSize',15);
  • lab7
% matched filter
% ‐ on‐off keying (OOK) using rectangular pulse (T=1)
% That is, when OFF, assuming A=0, thus E1 = 0
% when ON: assuming A=1, thus E2 = 1
% The "average bit energy" (Eb) = (E1+E2)/2 = 1/2
%
clear,clc,close all;
%% parameters
data_number = 10^5; % # of bits
EbN0dB_vec = [0:2:10 11.5]; % Eb/N0 in dB
Fs = 10; % sampling frequency (used to generate received samples)

%% transmitter
Data_bit = (rand(1,data_number) > 0.5 ); % random bits
p1 = ones(1,Fs); % discrete‐time rectangular pulse that represents one symbol
Data_pulse_array = (p1')*Data_bit;
Data_pulse = reshape(Data_pulse_array,1,length(p1)*data_number);

%% AWGN channel and receiver
% AWGN channel
[a, b] = size(Data_pulse);
% EbN0dB = 3;
EbN0dB = 10;
EbN0 = 10^(EbN0dB/10); % EbN0 is now in linear scale
sgma = sqrt(0.5/EbN0/2*10);
noise = normrnd(0, sgma, a, b);

Data_receive = Data_pulse + noise; % received samples
% receiver
D_filtered = conv(Data_receive,p1); % MF output
D_demapping = D_filtered(10:10:end) / 10; % sampling at symbol rate
% decsion based on D_demapping
D_demap_N = (D_demapping > 0.5); % >0.5: 1; <=0.5: 0

%% split two parts
index_0 = 1;
for i = 1:data_number
    if Data_bit(i) == 0 % TX == 0
        D_demapping_0(index_0) = D_demapping(i);
        index_0 = index_0 + 1;
    end
end

index_1 = 1;
for i = 1:data_number
    if Data_bit(i) == 1
        D_demapping_1(index_1) = D_demapping(i);
        index_1 = index_1 + 1;
    end
end

%% generate plots
figure;
hold on;
xaxis = -2.5:0.01:2.5;
histogram(D_demapping_0, 'FaceColor', '#0072BD','Normalization','pdf');
histogram(D_demapping_1, 'FaceColor', '#D95319','Normalization','pdf');

plot(xaxis,normpdf(xaxis,0,sgma / sqrt(10)), 'LineWidth',3); %need to scaling 10 with line 30
plot(xaxis,normpdf(xaxis,1,sgma / sqrt(10)),'r', 'LineWidth',3); %need to scaling 10 with line 30

xlabel('y(T)');
ylabel('probability density');
legend('bit 0 (simulation)', 'bit 1 (simulation)', 'bit 0 (theory)','bit 1 (theory)','FontSize',13);
grid on;
title('conditional PDF of (y(T) | bit 0 sent) and f(y(T) | bit 1 sent)','FontSize',15);