信号选取#

设被采样信号连续时间函数为 f(t) = 2sin(2pif1t) + 10sin(2pif2t) + 5sin(2pif3t) 其中， f1 = 100;
f2 = 400;
f3 = 800;

依据采样定理，设
fs = 2400;
Ts = 1/fs;

连续时间t和采样时间的关系可以表示为
t = n*Ts
其中n为序列的下标

取观察时间T = 0.1s;

使用matlab生成离散时间信号序列#

1
% 被采样信号设为f(t) = 2*sin(2*pi*f1*t) + 10*sin(2*pi*f2*t) + 5*sin(2*pi*f3*t)
2
f1 = 100;
3
f2 = 400;
4
f3 = 800;
5
fs = 2400;
6
Ts = 1/fs;
7
% 观察时间
8
T = 0.1;
9

10
% 序列索引
11
n = 0:round(T/Ts);
12

13
% 连续时间t和采样时刻的关系
14
t = n*Ts;
15

16
% 序列
17
x = 2*sin(2*pi*f1*t) + 10*sin(2*pi*f2*t) + 5*sin(2*pi*f3*t);
18
N = length(x);

标准DFT函数#

1
function X_w = DFT(x,N)
2
    for k = 0 : N - 1
3
        S = 0;
4
        for n = 0 : N - 1
5
            S = S + x(n+1) * exp(-1j*(2*pi*k*n/N));
6
        end
7
        X_w(k+1) = S;
8
        % X_clean = round(real(X_w), 12) + 1i*round(imag(X_w), 12);
9
    end
10
end

观察信号的频谱#

1
X_w = DFT(x, length(x));
2
k = 0:length(X_w)-1;
3

4
% 绘图
5

6
figure;
7

8
subplot(2,1,1);
9
stem(n, x, 'filled');
10
xlabel('n');
11
ylabel('x[n]');
12
title('时域序列');
13
grid on;
14

15
subplot(2,1,2);
16
stem(k, abs(X_w), 'filled');
17
xlabel('k');
18
ylabel('|X[k]|');
19
title('DFT幅度谱');
20
grid on;

可以得出：

现在的目标就是抑制频率为800hz的信号成分
可以设理想滤波器的通带截至频率为wp = 2 * pi * 500hz，阻带频率为ws = 2 * pi * 600hz

理想滤波器公式：

注意，公式中的角频率是不是Hz为单位，而是离散时间角频率，即rad/sample

连续时间信号 $x(t)=\cos(2\pi f t)$
采样后： $t=nT_s=\frac{n}{f_s}$
所以：

x[n] =\cos\left(2\pi f\frac{n}{f_s}\right) =\cos(\omega n)

因此不妨直接设
wp = wp = 2 * pi * 500 / fs
ws = 2 * pi * 600 / fs

则有截止频率
wc = (ws+wp) * 0.5;

如下matlab代码的理想滤波器实现

1
ws = 2*pi*500/fs;
2
wp = 2*pi*600/fs;
3
wc = (ws+wp) * 0.5;
4
alpha = (N - 1) / 2;
5
% 理想滤波器
6
% hd = sin(wc * (n - alpha)) ./ (pi * (n - alpha));
7
for m = 0 : N - 1
8
    if (n(m + 1) - alpha) ~= 0
9
        hd(m + 1) = sin(wc * (n(m + 1) - alpha)) / (pi * (n(m + 1) - alpha));
10
    else
11
        hd(m + 1) = wc / pi;
12
    end
13
end

(n(m + 1) - alpha) ~= 0表示公式中n - alpha不为0的时候，当n - alpha = 0时，可求得 hd = wc / pi

选择窗函数—Hann 窗为例#

1
% 窗函数
2
w = 0.5 * (1 - cos(2*pi*n / (N - 1)));

将理想滤波器与窗函数在时序相乘对理想滤波器进行截断，得到滤波器系数

1
% 滤波器
2
h = hd .* w;

查看滤波器的频域响应是否符合预期

1
Hw = DFT(h,length(h));
2
wk = 0 : length(h) - 1;
3
figure;
4

5
subplot(1,1,1);
6
stem(wk, abs(Hw), 'filled');
7
xlabel('wk');
8
ylabel('|hw|');
9
title('滤波器');
10
grid on;

由于DFT的周期性质，因此DTF图像上后部分有一段重复的图像。

得到滤波后的信号#

要得到滤波后的信号则直接将信号和滤波器进行卷积

1
x_f = conv(x, h); % 卷积之后长度会变化
2
n_xf = 0 : length(x_f) - 1;

注意卷积之后长度会变化，不再是N
再对滤波后的信号进行DTF分析

1
Xf_w = DTF(x_f, length(x_f));
2
k_xf = 0 : length(Xf_w) - 1;

最后绘制图像

1
figure;
2

3
subplot(3,1,1);
4
stem(wk, abs(Hw), 'filled');
5
xlabel('wk');
6
ylabel('|hw|');
7
title('滤波器');
8
grid on;
9

10
subplot(3,1,2);
11
stem(n_xf, x_f, 'filled');
12
xlabel('n_xf');
13
ylabel('x_f');
14
title('滤波后信号');
15
grid on;
16

17
subplot(3,1,3);
18
stem(k_xf, abs(Xf_w), 'filled');
19
xlabel('k_xf');
20
ylabel('|Xf_w|');
21
title('滤波后频谱');
22
grid on;

由于群延迟的存在，因此滤波后的时域信号产生了延迟，在图像上表现为往后延迟了，前100+个采样点的值在0附近。
由： $\phi(\omega)=-\frac{N-1}{2}\omega$
可得群延迟:
$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=\frac{N-1}{2}$
因此所有通过的频率分量大约都延迟：
$\frac{N-1}{2}$
个采样点，滤波后波形整体向右移动。

音乐

音乐

信号选取#

使用matlab生成离散时间信号序列#

标准DFT函数#

观察信号的频谱#

选择窗函数—Hann 窗为例#

得到滤波后的信号#

文章分享

音乐

窗函数法设计FIR滤波器

信号选取#

使用matlab生成离散时间信号序列#

标准DFT函数#

观察信号的频谱#

选择窗函数—Hann 窗为例#

得到滤波后的信号#

文章分享