在计算机飞速发展的今天,单纯用硬件进行信号接收与处理也变得效率低下,如何利用软件即各种优秀的算法对信号进行处理从而实现各种各样的功能,这已成为通信领域几乎必备的知识。今天我们来对一个某高校真实的例程进行剖析,进而更深刻的理解软件无线电。
先看任务内容:某软件无线电对0~5MHz信号进行采样,采样率为16MHz。已知在2MHz处有一带宽为40kHz的信号,设计信道化方案,要求最终输出的零中频信号采样率为100KHz。使用MATLAB对所设计的信道化方案进行仿真,并使用提供的采样数据对方案进行验证。
任务要求如下:
- 设计包括正交下变频和多级滤波抽取在内的信道化方案,要求使用CIC+HB+FIR滤波器的多级级联抽取,确定每一级的抽取倍数。
- 利用filterDesigner设计每一级HB/FIR滤波器,确定滤波器阶数,得到幅频响应和相频响应曲线,并导出滤波器系数。
- 编写m脚本文件,读取data.mat数据(16MHz采样率)作为系统输入,使用正交下变频、CIC(按所设计的阶数)以及第(2)步获得的各HB/FIR滤波器系数完成信道化仿真,以验证方案的正确性。
首先我们根据要求,零中频的信号采样率为100k,原信号采样率为16M,可知降采样倍数为16M/100k=160,而160分解因式可为5*2*2*2*2*2,由此可见,我们需要设计一个抽取比为5的CIC滤波器和5个半带滤波器,但是CIC一般不用做最后一级滤波,故最后用一个FIR滤波器收尾。
明确思路后,我们开始敲代码:
首先读取数据:
load data.mat; %% 读取原始采样数据
N = length(data); %% 数据长度
fs = 16e6; %% 原始采样率
Bw = 40e3; %% 待提取信号带宽
f_Lo = 2e6; %% 待提取信号中心频率
Power_xdBm(data(1:2^20),fs); %画原始采样信号功率谱图,data为从data.mat中读出的原始采样数据
title('输入信号功率谱图 采样率16MHz 频谱分辨率15.259Hz');
然后我们通过正交混频将频谱信号搬到零频处。
Lo = exp( - 1j * 2 * pi * f_Lo * ( 0 : N - 1 ) / fs ); %% 用于正交混频的本振信号
mixing_data = data .* Lo;
Power_xdBm_complex(mixing_data,fs); %画正交混频后信号功率谱图
title(['正交下变频信号功率谱图 采样率是',num2str(fs/1e6),'MHz 频谱分辨率15.259Hz']);
接下来进行CIC滤波器的设计,在设计滤波器时,可以用对应函数直接写,也可以在命令行输入filterDesigner进行窗口化设计。这里我们直接敲代码:
先设计其CIC滤波函数:
function output_data = cic_filter(N, fs, input_data, flag)
%%% CIC抽取滤波器
%%% 输入
% N:滤波器阶数,即抽取比
% fs:采样频率
% input_data:输入序列
% flag:画图标志位
%%% 输出
% output_data:滤波后的序列
b = zeros(1,N+1);b(1) = 1;b(N+1) = -1;
a = [1 -1];
[Hf,w] = freqz(b,a,[0:pi/1024:pi]);
output_data = filter(b,a,input_data);
end
将其运用到数据滤波中:
N = 5;
fs1 = fs/N;
cic_data = cic_filter(N,fs,mixing_data,0);
cic_data = cic_data(N+1:N:end);
Power_xdBm_complex(cic_data,fs1); %画CIC滤波后信号功率谱图
title(['CIC输出信号功率谱图 采样率是',num2str(fs1/10^6),'MHz 频谱分辨率15.259Hz']);
然后我们再设计第一级的HB滤波器:
function Hd = HBfilter1
%HBFILTER Returns a discrete-time filter object.
% FIR Window Halfband lowpass filter designed using the FIRHALFBAND
% function.
% All frequency values are in kHz.
Fs = 1600; % Sampling Frequency
Fpass = 50; % Passband Frequency
Dpass = 0.00057564620966; % Passband Ripple
% Calculate the coefficients using the FIR1 function.
b = firhalfband('minorder', Fpass/(Fs/2), Dpass, 'kaiser');
Hd = dfilt.dffir(b);
按照每次采样率下降一半的思路,我们可以设计出剩余四级滤波器。最后加上FIR滤波器:
function Hd = FIRfilter
%FIRFILTER Returns a discrete-time filter object.
% FIR Window Lowpass filter designed using the FIR1 function.
% All frequency values are in kHz.
Fs = 100; % Sampling Frequency
Fpass = 20; % Passband Frequency
Fstop = 25; % Stopband Frequency
Dpass = 0.057501127785; % Passband Ripple
Dstop = 0.0001; % Stopband Attenuation
flag = 'scale'; % Sampling Flag
% Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fpass Fstop]/(Fs/2), [1 0], [Dstop Dpass]);
% Calculate the coefficients using the FIR1 function.
b = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);
在设计完所有的滤波器后,我们将其加入到主函数中,便可实现对信号的160倍降采样处理,同时抽取出任务要求的2MHz附近带宽为40k的信号。
最后主函数如下:
clear ;
close all;
clc;
load data.mat; %% 读取原始采样数据
N = length(data); %% 数据长度
fs = 16e6; %% 原始采样率
Bw = 40e3; %% 待提取信号带宽
f_Lo = 2e6; %% 待提取信号中心频率
Power_xdBm(data(1:2^20),fs); %画原始采样信号功率谱图,data为从data.mat中读出的原始采样数据
title('输入信号功率谱图 采样率16MHz 频谱分辨率15.259Hz');
%% NCO mixing 正交混频
Lo = exp( - 1j * 2 * pi * f_Lo * ( 0 : N - 1 ) / fs ); %% 用于正交混频的本振信号
mixing_data = data .* Lo;
Power_xdBm_complex(mixing_data,fs); %画正交混频后信号功率谱图
title(['正交下变频信号功率谱图 采样率是',num2str(fs/1e6),'MHz 频谱分辨率15.259Hz']);
%% 第一级CIC滤波+抽取
N = 5;
fs1 = fs/N;
cic_data = cic_filter(N,fs,mixing_data,0);
cic_data = cic_data(N+1:N:end);
Power_xdBm_complex(cic_data,fs1); %画CIC滤波后信号功率谱图
title(['CIC输出信号功率谱图 采样率是',num2str(fs1/10^6),'MHz 频谱分辨率15.259Hz']);
%% HB1 Filter+2倍抽取
Hd1 = HBfilter1;
fs2 = fs1/2;
y1 = filter(Hd1,cic_data); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
y1 = y1(1:2:end);
Power_xdBm_complex(y1,fs2); %画HB1滤波后信号功率谱图
title(['HB1输出信号功率谱图 采样率是',num2str(fs2/1e6),'MHz 频谱分辨率15.259Hz']);
%% HB2 Filter+2倍抽取
Hd2 = HBfilter2;
fs3 = fs2/2;
y2 = filter(Hd2,y1); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
y2 = y2(1:2:end);
Power_xdBm_complex(y2,fs3); %画HB2滤波后信号功率谱图
title(['HB2输出信号功率谱图 采样率是',num2str(fs3/1e6),'MHz 频谱分辨率15.259Hz']);
%% HB3 Filter+2倍抽取
Hd3 = HBfilter3;
fs4 = fs3/2;
y3 = filter(Hd3,y2); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
y3 = y3(1:2:end);
Power_xdBm_complex(y3,fs4); %画HB3滤波后信号功率谱图
title(['HB3输出信号功率谱图 采样率是',num2str(fs4/1e6),'MHz 频谱分辨率15.259Hz']);
%% HB4 Filter+2倍抽取
Hd4 = HBfilter4;
fs5 = fs4/2;
y4 = filter(Hd4,y3); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
y4 = y4(1:2:end);
Power_xdBm_complex(y4,fs5); %画HB3滤波后信号功率谱图
title(['HB4输出信号功率谱图 采样率是',num2str(fs5),'Hz 频谱分辨率15.259Hz']);
%% HB5 Filter+2倍抽取
Hd5 = HBfilter5;
fs6 = fs5/2;
y5 = filter(Hd5,y4); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
y5 = y5(1:2:end);
Power_xdBm_complex(y5,fs6); %画HB3滤波后信号功率谱图
title(['HB5输出信号功率谱图 采样率是',num2str(fs6),'Hz 频谱分辨率15.259Hz']);
%% 最后一级FIR滤波+抽取
Hdfir = FIRfilter;
y6 = filter(Hdfir,y5); % 直接使用设计好的滤波器进行滤波,filter函数是滤波函数
Power_xdBm_complex(y6,fs6); %画HB3滤波后信号功率谱图
title(['FIR输出信号功率谱图 采样率是',num2str(fs6),'Hz 频谱分辨率15.259Hz']);
fs_out = fs6;
fir_data = y6;
%% 画出每一级滤波器的幅频响应
cic_data = cic_filter(N,fs,mixing_data,1);
title('CIC滤波器的幅频响应');
b1 = Hd1.Numerator;
figure;
freqz(b1,1,1024,fs2);
title('HB1滤波器的幅频响应');
b2 = Hd2.Numerator;
figure;
freqz(b2,1,1024,fs3);
title('HB2滤波器的幅频响应');
b3 = Hd3.Numerator;
figure;
freqz(b3,1,1024,fs4);
title('HB3滤波器的幅频响应');
b4 = Hd4.Numerator;
figure;
freqz(b4,1,1024,fs5);
title('HB4滤波器的幅频响应');
b5 = Hd5.Numerator;
figure;
freqz(b5,1,1024,fs6);
title('HB5滤波器的幅频响应');
bfir = Hdfir.Numerator;
figure;
freqz(bfir,1,1024,fs6);
title('FIR滤波器的幅频响应');
%% 画输出信号功率谱,其中fir_data为最后一级滤波抽取输出数据,fs_out为最后一级输出数据采样率
Power_xdBm_complex(fir_data(1:4096),fs_out); %
title('DDC输出信号功率谱图 采样率100kHz 频谱分辨率24.414Hz');
我们运行程序,可看到原信号频谱:
混频后频谱:
CIC滤波后频谱:
HB1滤波后频谱:
之后就是每级HB进行滤波后的频谱,我们不再赘述,直接看最后效果:
至此,我们便完成了任务要求,(伸个懒腰,哈~)。好了,如果你觉得这篇文章有用或者想要原始数据和代码及相关指导的,可以点击下方链接给鄙人打赏一手,我看到了都会回复的,嘻嘻,最后大家一起加油,在变强的路上一去不复返!
https://i.postimg.cc/nVRsQvvh/2-1.jpg