我已经搜索了(在google和dsp.stackexchange上都找到了),但发现有矛盾的答案。我一直在玩信号,这是结果。我对此不太了解。这是具有每四秒钟一次采样频率的信号。我设计了一个从0.8 mHz到1 mHz的过渡带的数字低通滤波器,并对信号进行了滤波。然后,我还设计了一个具有相同过渡带的高通滤波器,并对信号进行了滤波。这是结果。
第一张图片显示原始信号为黑色,低通信号为蓝色。它们几乎彼此叠在一起,但并不完全相同。红色曲线是信号减去位于信号顶部的高通信号。
第二张图像只是放大的第一张图像,显示了什么发生。在这里,我们可以清楚地看到两者并不相同。我的问题是为什么?是关于我如何实现这两个过滤器的事情,还是与我的实现无关的理论上的事情?我对滤波器的设计并不了解很多,但我确实知道这是违反直觉的。这是重现所有这些的完整MATLAB代码。我正在使用filtfilt命令消除相位延迟。但是,在此需要指出的另一点是,过滤器未标准化。当我做sum(Hd.Numerator)时,低通得到0.9930,高通得到0.007。我看不出如何解决这个问题。是否应该以某种方式缩放输出,因为系数不等于1?这种缩放可能与此有关吗?
close all
clear all
clc
data = dlmread('data.txt');
Fs = 0.25; % Sampling Frequency
N = 2674; % Order
Fpass = 0.8/1000; % Passband Frequency
Fstop = 1/1000; % Stopband Frequency
Wpass = 1; % Passband Weight
Wstop = 1; % Stopband Weight
dens = 20; % Density Factor
% Calculate the coefficients using the FIRPM function.
b = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datalowpassed = filtfilt(Hd.Numerator,1,data);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs = 0.25; % Sampling Frequency
N = 2674; % Order
Fstop = 0.8/1000; % Stopband Frequency
Fpass = 1/1000; % Passband Frequency
Wstop = 1; % Stopband Weight
Wpass = 1; % Passband Weight
dens = 20; % Density Factor
% Calculate the coefficients using the FIRPM function.
b = firpm(N, [0 Fstop Fpass Fs/2]/(Fs/2), [0 0 1 1], [Wstop Wpass], {dens});
Hd = dsp.FIRFilter('Numerator', b);
sum(Hd.Numerator)
datahighpassed = filtfilt(Hd.Numerator,1,data);
figure
subplot(2,1,1)
plot(data,'-ko')
hold on
plot(datalowpassed,'-bo')
plot(data-datahighpassed,'-ro')
legend('Original Signal','Low-Passed','Signal - High-Passed')
subplot(2,1,2)
plot(data-datalowpassed,'-bo')
hold on
plot(datahighpassed,'-ro')
legend('Signal - Low-Passed','High-Passed')
#1 楼
通常,您不能简单地从原始信号中减去信号的低通滤波版本以获得高通滤波信号。原因如下。您实际上正在做的是实现具有频率响应的系统$$ H(\ omega)= 1-H_ {LP}(\ omega)\ tag {1} $$
其中$ H_ {LP}(\ omega)$是低通滤波器的频率响应。请注意,$ H_ {LP}(\ omega)$是一个复杂的函数。您可能想要的是
$$ | H(\ omega)| = \ big | 1- | H_ {LP}(\ omega)| \ big | \ tag {2} $$
,但是通常情况下,当满足(1)时不是这种情况。
现在将$ H_ {LP}(\ omega)$写为
$ $ H_ {LP}(\ omega)= | H_ {LP}(\ omega)| e ^ {j \ phi(\ omega)} $$
其中$ \ phi(\ omega)$是低通滤波器的相位响应。如果您将过滤器定义为
$$ H_ {HP}(\ omega)= e ^ {j \ phi(\ omega)}-H_ {LP}(\ omega)= e ^ {j \ phi(\ omega)} \ left(1- | H_ {LP}(\ omega)| \ right)\ tag {3} $$
,那么这个新过滤器确实可以满足幅度关系(2)。因此,如果按照(3)进行减法,则可以实现幅度响应与低通滤波器的响应互补的高通滤波器。这意味着您需要先通过具有频率响应$ e ^ {j \ phi(\ omega)} $的全通滤波器对信号进行滤波,然后再从中减去信号的低通滤波版本。
实际上,如果低通滤波器具有线性相位响应,则这非常简单,因为相位项由
$$ e ^ {j \ phi(\ omega)}给出= e ^ {-j \ omega \ tau} \ tag {4} $$
,对应于一个简单的延迟。 FIR滤波器可以具有精确的线性相位响应,而由Parks-McClellan算法设计的滤波器确实具有线性相位。将(4)中的延迟$ \ tau $设为整数是有利的,这可以通过选择偶数滤波器阶$ n $来实现:$ \ tau = n / 2 $
必须要做的是:
设计偶数阶的线性相位FIR低通滤波器
使用
filter()
过滤信号;我们称结果$ x_ {LP} [n] $ 使原始信号延迟$ \ tau = n / 2 $个样本;我们称之为延迟信号$ x_d [n] $
通过减法生成高通滤波信号:$ x_ {HP} [n] = x_d [n] -x_ {LP} [n] $
这是Matlab / Octave中的一个非常简单的示例
h_lp = fir1(100,.3); % low-pass design h_hp = [zeros(50,1);1;zeros(50,1)] - h_lp; % high-pass design by subtraction [H_lp,w] = freqz(h_lp,1,1024); [H_hp,w] = freqz(h_hp,1,1024); plot(w/2/pi,20*log10(abs(H_lp)),w/2/pi,20*log10(abs(H_hp))) grid, axis([0,.5,-100,5])
编辑:
由于使用了
filtfilt
命令,因此人为消除了相位,并且上面的等式(1)和(2)变得等效,因为所有频率响应实际上都是设计响应的平方幅度。因此,除了由于滤波器设计过程而产生的细微差异,数值误差以及由filtfilt
函数引起的细微差异(自动选择的初始条件使开始和结束瞬态最小化)外,从原始数据中减去滤波后的数据的结果应与直接相似。用互补滤波器进行滤波。由于在您的示例中情况并非如此,因此我怀疑由于极高的滤波器阶数,滤波器设计例程会给您带来麻烦。如果使用更简单的过滤器(我选择了$ n = 100 $)执行相同的操作,则可以得到期望的结果。在下图中,您会看到蓝色的数据部分,绿色的低通滤波器的输出,以及红色的原始数据减去高通滤波器的输出的结果。绿色和红色曲线几乎完全相同。x = load('data.txt'); % data to be filtered h_lp = fir1(100,.3); % LP impulse response h_hp = fir1(100,.3,'high'); % HP impulse response y = filtfilt(h_lp,1,x); % apply low pass filter yh = filtfilt(h_hp,1,x); % apply high pass filter yd = x - yh; % low pass by difference with high pass filter n = 1:length(x); plot(n,x,n,y,'g.',n,yd,'r') axis([3500,4000,140,150])
评论
$ \ begingroup $
如果您试图以此方式设计高通滤波器,则必须注意低通滤波器的规格。低通滤波器中的阻带衰减通常足够高,可以在高通中实现较小的通带纹波,但是LP滤波器中的通带纹波通常无法在HP滤波器中实现足够的阻带衰减。
$ \ endgroup $
–大卫
15年3月6日在13:49
$ \ begingroup $
感谢您的详细回复。它清除了几件事。
$ \ endgroup $
–固定点
2015年3月7日,0:55
#2 楼
关于缩放比例:如果对系数求和,将得到DC的幅度。因此,获得这些数字是有意义的(LP为$ \大约1 $,HP为$ \大约0 $)。
除了Matt L.的出色答案外,您还可以指出他使用的被称为幅度互补滤波器,这是线性相位FIR滤波器的常见情况,即,
$ | H_ {LP} | + | H_ {HP} | = 1 $
当从两个并行的全通部分创建滤波器并增加/减去输出时,低通/高通滤波器将是功率互补的,即,
$ | H_ {LP} | ^ 2 + | H_ {HP} | ^ 2 = 1 $
评论
您的过滤器订单非常高。您可能会在设计过程中遇到数字问题。通过计算和绘制FFT的幅度来检查设计的滤波器。此外,请参阅下面的答案,了解如何使用减法从低通滤波器生成高通滤波器。您可以通过这种方式制作高通信号,但滚降的顺序始终是一阶:sound.westhost.com/articles/derived-xovers.htm如果您先将信号延迟LPF的群延迟,然后再减去,您可以得到三阶