Matlab里關于T2F/F2T/lpf.m函式呼叫應注意的問題及函式修改

最近看到網上很多人在做基本的數字信號調制解調仿真時候用到了T2F/F2T/lpf.m這三個流傳甚廣的函式，這三個函式在信號解調部分組成了一個抽樣判決前的低通濾波器，以前我也直接呼叫過，但是今天打開函式一看，還是發現了點問題，先看看T2F函式：

1. T2F.m

function [f,sf]= T2F(t,st)
%Input is the time and the signal vectors,the length of time must greater
%than 2
%Output is the frequency and the signal spectrum
dt = t(2)-t(1);
T=t(end);
df = 1/T;
N = length(st);
f=-N/2*df : df : N/2*df-df;
sf = fft(st);
sf = T/N*fftshift(sf);

T2F.m其實就是計算輸入序列的FFT，輸入引數為時間序列t和信號序列st，回傳值為頻率序列f和對應的幅值sf，函式先計算了時間采樣間隔dt，然后將時間序列t的最后一位（截至時刻）賦值給了T，接著定義頻率采樣間隔df=1/T，采樣點數N，然后進行N點的fft變換→定義頻率序列→呼叫fftshift，

正常來說，大部分人撰寫信號，時間序列的起始為0時刻，這樣的話T=t(end)=(N-1)dt，這種情況下該函式在定義頻率采樣間隔時：df=1/T=1/(N-1)/dt.很明顯該結果是錯誤的，因為采樣頻率fs=1/dt，df=fs/N=1/N/dt，

這樣會使得當設計的信號采樣點數減小的情況下，計算的頻譜產生很大的偏差，（但是如果輸入信號時間序列的起始為dt時刻的話，呼叫該函式是正確的，但這種情況應該比較少吧）

將該函式稍作修改就能用啦，如下（fftshift前面的系數沒什么用我就刪了）：

function [f,sf]= T2F1(t,st)
%需滿足時刻起始點為0
%Input is the time and the signal vectors,the length of time must greater
%than 2
%Output is the frequency and the signal spectrum
dt = t(2)-t(1);
T=t(end);
df = 1/(T+dt);
N = length(st);
f=-N/2*df : df : N/2*df-df;
sf = fft(st);
sf = fftshift(sf);

2. lpf.m

function [t,st]=lpf(f,sf,B)
%This function filter an input data using a lowpass filter
%Inputs: f: frequency samples
% sf: input data spectrum samples
% B: lowpass bandwidth with a rectangle lowpass
%Outputs: t: time samples
% st: output data time samples
df = f(2)-f(1);
T = 1/df;
hf = zeros(1,length(f));
bf = [-floor( B/df ): floor( B/df )] + floor( length(f)/2 );
hf(bf)=1;
yf=hf.*sf;
[t,st]=F2T(f,yf);
st = real(st);

接著看lpf.m函式，呼叫這個函式要輸入引數為頻率序列f，對應的幅值sf，矩形低通濾波器截至頻率B，采樣頻率間隔df，然后計算T（顯然也有問題，但是該賦值后續并沒有使用，因此不管它），hf為長度與f相同的全零序列，而bf顯然是想計算f中對應低通頻段的位置，我們知道f中第位表示的是直流（f=0），那么bf的運算式里少加了一個1，而且floor是向下取整，我們應該改為fix，即向0取整，最后一行取實部也應該不必要，如果原始的時間序列是實的，那么st也是實序列，

修改后的函式如下：

function [t,st]=lpf1(f,sf,B)
%This function filter an input data using a lowpass filter
%Inputs: f: frequency samples
% sf: input data spectrum samples
% B: lowpass bandwidth with a rectangle lowpass
%Outputs: t: time samples
% st: output data time samples
df = f(2)-f(1);
T = 1/df;
hf = zeros(1,length(f));
bf = [-fix( B/df ): fix( B/df )] + floor( length(f)/2 )+1;
hf(bf)=1;
yf=hf.*sf;
[t,st]=F2T1(f,yf);
st = real(st);

3.F2T

function [t,st]=F2T(f,sf)
%This function calculate the time signal using ifft function for the input
df = f(2)-f(1);
Fmx = ( f(end)-f(1) +df);
dt = 1/Fmx;
N = length(sf);
T = dt*N;
t = 0:dt:T-dt;
sff = fftshift(sf);
st = Fmx*ifft(sff);

F2T.m函式就是計算信號的IFFT，基本沒什么問題，Fmx算的就是采樣頻率，寫得這么奇怪應該是別人修改后的版本，而且時間序列的運算式在這又是明顯正確的，也稍微改一下，一樣把ifft前面的系數刪去，如下：

function [t,st]=F2T1(f,sf)
%This function calculate the time signal using ifft function for the input
df = f(2)-f(1);
N = length(sf);
Fmx = N*df;%采樣頻率
dt = 1/Fmx;%采樣間隔
T = dt*N;
t = 0:dt:T-dt;
sff = fftshift(sf);
st = ifft(sff);

原始的三個函式在信號序列采樣點數較大時對解調影響較小，親測PSK和FSK能正確解調，但是當采樣點數只有幾十個時，解調誤碼率很高，因為組成的低通濾波器的頻帶完全錯位和擴寬了，以上修改后的函式親測能運行且正確解調，

最后再補充一句，呼叫現成的函式時一定一定一定要注意適用條件，能瞅一眼就都瞅一眼，以免后續產生各種問題，