在 matlab 上设计希尔伯特滤波器(纯 90 相位和幅度 0db)

信息处理 matlab 过滤器设计 阶段 希尔伯特变换
2022-02-07 05:41:32

我需要设计一个能够测量交流电并将输入相位偏移 90° 的滤波器。

滤波器应该在 40 到 60Hz 之间工作,并且在这个范围内应该有单位增益。在这个范围之外,我不关心增益或相位。

我研究了希尔伯特滤波器,听起来这种滤波器是我需要的,但我不知道如何使用。

让我粘贴我的 matlab 代码(在马特推荐之后修复):

%     close all
clear all
clc

%Period
Fs = 1000; %sample frequency

fl = 40;   %lower frequency 

f_min=fl/(Fs/2);

b = firpm(40,[f_min (1-f_min)],[1 1],'h');   % Bandpass Hilbert
%     fvtool(b)

t = 0:1/Fs:1.2;
t_window=0.2;

%input signal Frequency

f=40;
%Generate AC signal (input)- at 40hz (lower limit)
y1 = sin(2*pi*f*t);

%Generate filtered signal-40hz
y1f=filter(b,1,y1);
y1 = [zeros(1,(length(b)-1)/2), y1(1:end-(length(b)-1)/2)];

f=60;
%Generate AC signal (input)- at 60Hz (upper limit)
y2 = sin(2*pi*f*t);

%Generate filtered signal-60hz
y2f=filter(b,1,y2);
y2 = [zeros(1,(length(b)-1)/2), y2(1:end-(length(b)-1)/2)];

figure
hold on

plot(t,y1,'*-')
plot(t,y1f)
plot(t,y2,'*-')
plot(t,y2f)

hold off
grid

legend('input','filtered')
axis([max(t)-t_window max(t) -1.1 1.1])

range(y1f(end-(length(t)/10):end))/range(y2f(end-(length(t)/10):end))

我原以为我的信号yf会发生 90° 偏移,但这不起作用。

有人对如何解决这个问题以及我的代码有什么问题有一些想法?


第 2 部分 - 感谢 Hilmar 和 Richard。让我修复我的代码并使用您的代码并展示正在发生的事情。完毕。但这还行不通...

  • 目标 1- 90 相移 - 好!(多谢你们!)
  • 目标2-没有变化的幅度-不好。如果我将输入频率从 60 更改为 40Hz,信号会衰减 5%。

你能帮我解释一下为什么会发生这种情况吗?我应该怎么做才能解决这个问题?....我需要在这个频率之间尽可能平坦的幅度。我试图增加过滤器的顺序(40-> 60),但还没有解决我的问题。


第 3 部分 - 谢谢马特。让我修复我的代码并使用您的建议。现在工作!

现在,我的滤波信号得到了适当的偏移,并且在 40 到 60hz 之间具有低衰减。在此范围内衰减小于 0.4%。

再次感谢大家(马特、理查德和希尔玛)

4个回答

Hilbert Transformers are non-causal, i.e. they need to be delayed to be implementable. So you get the 90 degree phase-shift plus a bulk delay of 20 samples (half the filter length).

You see the 90 degree phase shift if you delay the original signal by 20 samples as well.

EDIT for Part 2:

Your lower bandpass cutoff is too high. It currently sits at 50 Hz. It needs to go down to 40 Hz or potentially a little lower. You may also have to crank up the number of points to get less amplitude ripple in the pass band. Something like

b = firpm(60,[38 950]*2/2000,[1 1],'h');

gets you about 0.1 dB ripple. Cranking the number of tabs up to 96 will reduce the pass band ripple to 0.01dB.

@Rodrigo PG:您的 Matlab Hilbert 变换器是因果关系。这意味着,为了及时同步您的y与你的顺序yf序列,你需要延迟你的y序列(长度为b减去 1)/2 个样本时的长度b很奇怪。这会延迟您的输入y by the same amount as the Hilbert transformer's inherent time delay. Insert the following line of code after your yf = filter(b,1,y); command:

y = [zeros(1,(length(b)-1)/2), y];

to see the desired 90-deg phase shift.

Regarding the magnitude response of the FIR Hilbert transformer, it is never ideally flat. The passband magnitude ripple can be reduced by increasing the filter order. However, in your case the frequency 40 (whatever units) at a sampling frequency of 2000 is outside the passband of the designed filter, because your lower band edge is 0.05 (which corresponds to a frequency of 50 for the given sampling rate).

So, if you want the Hilbert transformer to work well at a frequency of 40, this frequency must be inside the passband, i.e. you must lower the lower band edge to say 0.035. If the resulting passband ripple is too high, then you also must increase the filter order.

@Rodrigo PG: Now that your Hilbert filter is working, you can start experimenting.

[1] Window your Hilb filter’s coefficients with a Hanning (Hann) or Hamming window to reduce the passband ripple. But make sure your signal of interest is still within the Hilb filter’s passband.

[2] Decimate your input signal by a factor of two using a halfband filter. That will allow you to use a simpler (fewer coefficients) Hilb filter.

[3] Insert a zero-valued coefficient in between each of your original Hilb filter’s coefficients. That will enable you to perform a Hilbert transform on lower-frequency signals.