Matlab里关于T2F/F2T/lpf.m函数调用应注意的问题及函数修改

本文针对网络上广泛使用的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中第\frac{N}{2}+1位表示的是直流(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能正确解调,但是当采样点数只有几十个时,解调误码率很高,因为组成的低通滤波器的频带完全错位和扩宽了,以上修改后的函数亲测能运行且正确解调。

最后再补充一句,调用现成的函数时一定一定一定要注意适用条件,能瞅一眼就都瞅一眼,以免后续产生各种问题。

``` % 参数设置 fm = 5; % 基带信号频率 (Hz) fc = 100; % 载波频率 (Hz) A = 2; % 载波幅度 N0 = 0.01; % 噪声单边功率谱密度 B = 2 * fm; % 噪声带宽 fs = 10 * fc; % 采样频率 t = 0:1/fs:1-1/fs; % 时间向量 (1秒) % 信源信号 mt = sqrt(2) * cos(2 * pi * fm * t); % 基带信号 % AM调制 s_am = (A + mt) .* cos(2 * pi * fc * t); % AM调制信号 % 窄带高斯噪声生成 noise = noise_nb(fc, B, N0, t); % 生成窄带噪声 s_am_noisy = s_am + noise; % 添加噪声 % AM解调 rt = s_am_noisy .* cos(2 * pi * fc * t); % 相干解调 rt = rt - mean(rt); % 去除直流分量 % 傅叶变换 [f, rf] = T2F(t, rt); % 将解调信号转换到频域 % 低通滤波 [t, rt_filtered] = lpf(f, rf, fm); % 低通滤波,截止频率为 fm % 绘图 figure; subplot(3, 1, 1); plot(t, mt); title('基带信号 m(t)'); xlabel('时间 (s)'); ylabel('幅度'); subplot(3, 1, 2); plot(t, s_am_noisy); title('AM调制信号 + 噪声'); xlabel('时间 (s)'); ylabel('幅度'); subplot(3, 1, 3); plot(t, rt_filtered); title('解调后的基带信号'); xlabel('时间 (s)'); ylabel('幅度'); % 辅助函数:傅叶变换 function [f, sf] = T2F(t, st) dt = t(2) - t(1); T = t(end); df = 1 / T; N = length(st); f = (-N/2 : N/2-1) * df; % 频率向量 sf = fft(st); sf = fftshift(sf) * dt; % 傅叶变换结果 end % 辅助函数:低通滤波 function [t, st] = lpf(f, sf, B) df = f(2) - f(1); hf = zeros(1, length(f)); % 初始化滤波器 hf(abs(f) <= B) = 1; % 矩形低通滤波器 yf = hf .* sf; % 频域滤波 [t, st] = F2T(f, yf); % 逆傅叶变换 st = real(st); % 取实部 end % 辅助函数:逆傅叶变换 function [t, st] = F2T(f, sf) df = f(2) - f(1); N = length(sf); t = (0 : N-1) / (N * df); % 时间向量 st = ifft(ifftshift(sf)) * N * df; % 逆傅叶变换 end % 辅助函数:窄带高斯噪声生成 function [out] = noise_nb(fc, B, N0, t) dt = t(2) - t(1); Fmx = 1 / dt; % 采样频率 n_len = length(t); p = N0 * Fmx / 2; % 噪声方差 rn = sqrt(p) * randn(1, n_len); % 生成白噪声 [f, rf] = T2F(t, rn); % 转换到频域 [t, out] = bpf(f, rf, fc - B/2, fc + B/2); % 带通滤波 end % 辅助函数:带通滤波 function [t, st] = bpf(f, sf, f1, f2) df = f(2) - f(1); hf = zeros(1, length(f)); % 初始化滤波器 hf((f >= f1) & (f <= f2)) = 1; % 矩形带通滤波器 hf((f <= -f1) & (f >= -f2)) = 1; % 对称负频率部分 yf = hf .* sf; % 频域滤波 [t, st] = F2T(f, yf); % 逆傅叶变换 st = real(st); % 取实部 end```为什么这个matlab程序每次运行后解调后的基带图像都不一样
03-26
%% 基础 clear all close all i=10; %%基带信号码元数 j=5000; a=randi([0,1],1,i); %%产生随机序列 t=linspace(0,5,j); %0-5之间产生5000个点行矢量,即将[0,5]分成5000份(采样速率为1000) f1=5; f2=10; f3=15; f4=20; f_4fsk=[5,15,20,25]; fm=i/5; EsN0=0.1; %信噪比 %% 产生基带信号 st1=t; for n=1:10 if a(n)<1 for m=j/i*(n-1)+1:j/i*n st1(m)=0; end else for m=j/i*(n-1)+1:j/i*n st1(m)=1; end end end figure(1); subplot(221); plot(t,st1); axis([0,5,-1,2]); title('基带信号st'); %% 4进制信号生成 result=reshape(a,2,5);% 将二进制信号重塑为2行m列的矩阵,每列是一组2位二进制数 quaternary_signal=2*result(1,:)+result(2,:);% 转换公式:第1位(高位)×2 + 第2位(低位),得到0~3的四进制符号 num_symbols = length(quaternary_signal); samples_per_symbol = j / num_symbols; % 每个符号采样点数=1000 s_4fsk = zeros(1, j); % 总长度j=5000 for idx = 1:num_symbols % 当前符号的时间索引 start_idx = round((idx-1)*samples_per_symbol + 1); end_idx = round(idx * samples_per_symbol); % 选择频率 freq = f_4fsk(quaternary_signal(idx)+1); % 生成当前符号对的4FSK段 s_4fsk(start_idx:end_idx) = cos(2*pi*freq*t(start_idx:end_idx)); end subplot(222); plot(t,s_4fsk); axis([0,5,-2,2]); title('4FSK信号'); %% 加入噪声 fsk=awgn(s_4fsk,EsN0,'measured'); %高斯白噪声信道 subplot(223); plot(t,fsk); axis([0,5,-2,2]); title('加入噪声的信号'); %% 生成4个载波信号(保持原代码不变) carrier = cell(1, 4); for k = 1:4 carrier{k} = cos(2*pi*f_4fsk(k)*t); end %% 相干解调(完全保留原代码结构) demod = zeros(4, j); % 4路解调结果 for k = 1:4 mult = fsk .* carrier{k}; % 与对载波相乘 [freq, sf] = T2F(t, mult); % 傅叶变换 [t, demod(k, :)] = lpf(freq, sf, 2*fm); % 低通滤波 demod(k, :) = real(demod(k, :)); end % 采样判决(保留原循环和索引逻辑) at = t; len = 500; % 每个码元的采样点数(不变) mid = 250; % 码元中点偏移量(不变) for m = 0:9 mid_idx = m * len + mid; % 码元中点索引(不变) % 提取4路解调信号中点值(逻辑不变) val1 = demod(1, mid_idx); % 对5Hz载波 val2 = demod(2, mid_idx); % 对10Hz载波 val3 = demod(3, mid_idx); % 对15Hz载波 val4 = demod(4, mid_idx); % 对20Hz载波 % 判决四进制状态(逻辑不变) [~, state] = max([val1, val2, val3, val4]); state_1 = state - 1; % 转换为0-3 % 赋值当前码元区间(循环结构不变) % 判决最大相关值 for idx = m*len + 1 : (m+1)*len at(1, idx) = state_1; end end %% 抽样判决 subplot(224); plot(t,at); axis([0,5,-1,3]); title('抽样判决后波形') function [f,sf]= T2F(t,st) %利用FFT计算信号的频谱并与信号的真实频谱的抽样比较。 %脚本文件T2F.m定义了函数T2F,计算信号的傅立叶变换。 %Input is the time and the signal vectors,the length of time must greater %than 2 %Output is the frequency and the signal spectrum T=t(end); df = 1/T; %df实际上等于采样速率/总点数,即为频谱横轴步进值 N = length(st); f=-N/2*df : df : N/2*df-df; %跟奈奎斯特采样定理有关,f取值范围在-fs到fs之间,fs为采样速率(信号最大频率不超过采样速率的一半) % f=(-N/2 : N/2-1)/T; sf = fft(st); sf = T/N*fftshift(sf); end function [t,st]=F2T(f,sf) %脚本文件F2T.m定义了函数F2T,计算信号的反傅立叶变换。 %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=-T/2:dt:T/2-dt; t = 0:dt:T-dt; sff = fftshift(sf); st = Fmx*ifft(sff); end 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); end 对该程序生成流程图
最新发布
11-28
评论 11
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值