学习随笔之数字陷波器

进行仿真实验，从某多频率信号中提取指定频率信号，并计算其幅值。具体如下：

生成仿真信号x0，含3频率f0、f1、f2，其中f0为待提取频率；
进行陷波器设计，在指定频率f0与采样频率fs计算陷波数字频率；
将原始信号x0通过陷波器，去除f0频率成分，生成中间信号x1；
将原始信号x0减去中间信号x1，得到只包含f0成分的最终信号x2；
将x2的一个采样周期(f0)进行面积积分，利用正弦波面积积分法计算其幅值（对于正弦波X=Asin(2πft)，其单个整周期面积为2A/fπ）。

仿真代码：

%% 原始信号
fs = 500;       % 采样频率
Ts = 1/fs;       % 采样周期
t = Ts:Ts:600;   % 采集10min的数据
len = length(t);  % 数据长度
f0 = 1;  A0 = 10;% 三个频率成分，f0为要提取的成分,A0为对应幅值待测
f1 = 5;  A1 = 2;
f2 = 10; A2 = 1;
x0 = A0*sin(2*pi*f0*t)+A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t);  % 原始信号（由三个频率成分组成,f0为待提取成分）
%% 陷波器设计
w0 = 2*pi*f0/fs;  % 陷波器数字频率（弧度单位）
r = 0.999;        % 滤波常数，决定陷波器阻带与衰减特性
b_iir = [1 -2*cos(w0) 1];
a_iir = [1 -2*r*cos(w0) r*r];
figure            % 陷波器幅频特性与相频特性
freqz(b_iir,a_iir,1024*16)
%% 陷波器滤波
x1 = filter(b_iir,a_iir,x0); %陷波器滤波，去除f0成分
x2 = x0 - x1;      % 滤波结果同原信号相减，提取f0成分
%% 滤波结果分析
NFFT = 2^(nextpow2(len)+5); % 进行FFT变换所用数，此处的+5为提高频谱分辨率所加
Y0 = fft(x0,NFFT)/len;      % 进行FFT变换与归一化（滤波前信号）
Y2 = fft(x2,NFFT)/len;      % 进行FFT变换与归一化（滤波后信号）
f = fs/2*linspace(0,1,NFFT/2+1); % FFT结果对应频率值
% 时域
n1 = 1:10*fs;  % 显示前10s的时域数据
figure
plot(t(n1),x0(n1),t(n1),x2(n1))
title('\fontname{微软雅黑}\fontsize{10}时域')
xlabel('\fontname{微软雅黑}\fontsize{10}时间(s)')
ylabel('\fontname{微软雅黑}\fontsize{10}波形')
legend('\fontname{微软雅黑}\fontsize{10}滤波前','\fontname{微软雅黑}\fontsize{10}滤波后')
grid on
% 频域
n2 = find(f<15); % 显示0-15Hz的频域数据
figure
plot(f(n2),2*abs(Y0(n2)),'-o',f(n2),2*abs(Y2(n2)),'-*')
title('\fontname{微软雅黑}\fontsize{10}频域')
xlabel('\fontname{微软雅黑}\fontsize{10}频率(Hz)')
ylabel('\fontname{微软雅黑}\fontsize{10}幅值')
legend('\fontname{微软雅黑}\fontsize{10}滤波前','\fontname{微软雅黑}\fontsize{10}滤波后')
grid on
%% 提取f0成分幅值（面积积分法）
jj = 0;A = 0;sumx = 0;t2 = 0;
for ii = 1:len
    sumx = sumx + abs(x2(ii))*Ts; % 面积积分
    if rem(ii,1/f0*fs)==0  % 在每个整采样周期进行解算
        jj = jj+1;
        A(jj) = sumx*f0*pi/2; % 利用面积积分结果计算正弦波幅值
        t2(jj) = t(ii);       % 每次计算对应时刻
        sumx = 0;
    end    
end
figure
plot(t2,A);
title('\fontname{微软雅黑}\fontsize{10}面积积分法提取幅值')
xlabel('\fontname{微软雅黑}\fontsize{10}时间(s)')
ylabel('\fontname{微软雅黑}\fontsize{10}幅值')
grid on