Skip to content

LTI 系统频域分析

相位

我们考虑频率响应为 \(H\left(e^{j\omega \;} \right)=e^{-j\omega t_0 }\)

它的增益不必说, 是单位增益. 同时它有着线性相位, 即

\[ |H\left(j\omega \right)=1|,\angle H\left(j\omega \right)=\omega t_0 \]

我们当然知道, 这种频率响应产生的输出就是输入的时间移位, 即

\[ y\left(t\right)=x\left(t-t_0 \right) \]

对离散信号来说, 当移位的 \(n_0\) 是整数时, \(y\left\lbrack n\right\rbrack =x\left\lbrack n-n_0 \right\rbrack\) . 否则就没有这么简单了, 可见 dsp 中 problem4.7, 是包络的时间移位, 一个 sinc 函数.

有些系统虽然有着单位增益, 但是是非线性相位, 此时信号就会产生很大的变化. 当然, 对这种有单位增益的系统, 我们称之为全通系统. 它的幅频特性完全由相位特性决定.

群时延

前面说到, 对于线性相位系统来说, 相位特性的斜率就是时移的大小. 那么对于非线性相位系统来说, 我们也可以定义一个类似的群时延, 等于那个频率上相位特性斜率的负数, 表示不同频率的相位的时间移位.

\[ \tau \left(\omega \right)=-\frac{\mathrm{d}}{\mathrm{d}\omega }\left\lbrace \angle H\left(j\omega \right)\right\rbrace \]

example 6.1

考虑下面的频率响应 \(H\left(j\omega \right)\)

\[ H\left(j\omega \right)=\prod_{i=1}^3 H_i \left(j\omega \right) \]

其中

\[ H_{i\;} \left(j\omega \right)=\frac{1+{\left(j\frac{\omega }{\omega_i }\right)}^2 -2j\xi_i \left(\frac{\omega }{\omega_i }\right)}{1+{\left(j\frac{\omega }{\omega_i }\right)}^2 +2j\xi_i \left(\frac{\omega }{\omega_i }\right)} \]
\[ \left\lbrace \begin{array}{lll} \omega_1 =315\;{\mathrm{r}\mathrm{a}\mathrm{d}}/{\mathrm{s}}, & \xi_1 =0\ldotp 066\\ \omega_2 =943\;{\mathrm{r}\mathrm{a}\mathrm{d}}/{\mathrm{s}}, & \xi_2 =0\ldotp 033\\ \omega_3 =1888\;{\mathrm{r}\mathrm{a}\mathrm{d}}/{\mathrm{s}}, & \xi_3 =0\ldotp 058 \end{array}\right. \]

实际上上面的三个角速度分别对于 50 Hz, 150 Hz, 300 Hz 的频率.

它是一个全通系统, 我们不妨看一下它的相位特性. 

w1 = 315; w2 = 943; w3 = 1888;
xi1 = 0.066; xi2 = 0.033; xi3 = 0.058;
w = 0:0.01:3000;
H1 = (1+(1j*w/w1).^2-2*1j*xi1*(w/w1))./((1+(1j*w/w1).^2+2*1j*xi1*(w/w1)));
H2 = (1+(1j*w/w2).^2-2*1j*xi2*(w/w2))./((1+(1j*w/w2).^2+2*1j*xi2*(w/w2)));
H3 = (1+(1j*w/w3).^2-2*1j*xi3*(w/w3))./((1+(1j*w/w3).^2+2*1j*xi3*(w/w3)));
H = H1.*H2.*H3;
angle_H = angle(H);
plot(w./(2*pi), angle_H);
xlim([0,400]);
xlabel("frequency/Hz")
ylabel("phase(rad/s)")

figure_0.png

可以用 unwrap 函数平移 \(2k\pi\) 个相位使得相位连续:

u_angle_H = unwrap(angle_H);
plot(w./(2*pi), u_angle_H);
xlim([0,400]);
xlabel("frequency/Hz")
ylabel("phase(rad/s)")

figure_1.png

求群时延:

tau = [-diff(unwrap(angle_H)), 0];
plot(w./(2*pi), tau);
xlim([0,400]);
xlabel("frequency/Hz")
ylabel("group delay/s")

figure_2.png

输入中的不同频率被延时了不同的量, 这种现象称为弥散 (dispersion) . 这个例子中, 群时延在 50 Hz 处最大, 所以单位冲激响应的后面部分会在接近 50 Hz 的较低频率上振荡:

h = ifft(H);
% ps 后面基本上都是0, 没有必要画了, 只截取前面部分即可.
% ifft 变换回去后, 怎么确定时间? 
h1 = h(1:100);
plot(real(h1));

figure_3.png

对数模和相位图

在对数标尺上, 我们可以看到阻带的模特性细节.

伯德图: \(20\log_{10} |H\left(j\omega \right)|\) 和 \(\angle H\left(j\omega \right)\) 对于 \(\log_{10\;} \left(\omega \right)\) 的图称为 Bode 图.

matlab 可以很轻松地画出伯德图:

bode(sys1,sys2,...,sysN)

举个例子:

\[ H\left(j\omega \right)=\frac{100\left(1+j\omega \right)}{\left(10+j\omega \right)\left(100+j\omega \right)} \]

将之写成多个系统相乘的形式, 很容易证明最后的伯德图是各个伯德图的相加.

\[ H\left(j\omega \right)=\left(\frac{1}{10}\right)\left(\frac{1}{1+j\frac{\omega }{10}}\right)\left(\frac{1}{1+j\frac{\omega }{100}}\right)\left(1+j\omega \right) \]
H1 = tf([0.1,], [0.1, 1]);
H2 = tf([1,], [0.01, 1]);
H3 = tf([1,1], [1,]);
H = H1*H2*H3

H =

       0.1 s + 0.1
  ----------------------
  0.001 s^2 + 0.11 s + 1

连续时间传递函数。
模型属性

bode(H)

figure_4.png

二阶电路

再举一个二阶电路的例子:

\[ H\left(j\omega \right)=\frac{1}{{\left(j\frac{\omega }{\omega_n }\right)}^2 +2\xi \left(j\frac{\omega }{\omega_n }\right)+1} \]

下面取 \(\omega =200\pi \;{\mathrm{s}}^{-1} ,\xi =0\ldotp 1,0\ldotp 2,0\ldotp 4,0\ldotp 7,1,1\ldotp 5\) 做伯德图

w1 = 200;
xi1 = 0.1; xi2 = 0.2; xi3 = 0.4; xi4 = 0.7; xi5 = 1; xi6 = 1.5;
H1 = tf(1, [1,2*xi1,1]);
H2 = tf(1, [1,2*xi2,1]);
H3 = tf(1, [1,2*xi3,1]);
H4 = tf(1, [1,2*xi4,1]);
H5 = tf(1, [1,2*xi5,1]);
H6 = tf(1, [1,2*xi6,1]);

bode(H1,H2,H3,H4,H5,H6);
legend("\xi_1 = 0.1","\xi_2 = 0.2","\xi_3 = 0.4","\xi_4 = 0.7","\xi_5 = 1", "\xi_6 = 1.5")

figure_5.png

用 freqresp 函数可以得到 tf 函数的频率响应.

% 设置频率范围
frequencies = logspace(-3, 1, 1000);  % 生成对数均匀分布的频率点 10^-3到10^1

% 计算频率响应
resp1 = freqresp(H1, frequencies);
resp2 = freqresp(H2, frequencies);
resp3 = freqresp(H3, frequencies);
resp4 = freqresp(H4, frequencies);
resp5 = freqresp(H5, frequencies);
resp6 = freqresp(H6, frequencies);
% 绘制频率响应曲线
figure;
hold on;
%squeeze删除长度是1的维度
semilogx(frequencies, abs(squeeze(resp1))); 
semilogx(frequencies, abs(squeeze(resp2))); 
semilogx(frequencies, abs(squeeze(resp3))); 
semilogx(frequencies, abs(squeeze(resp4))); 
semilogx(frequencies, abs(squeeze(resp5))); 
semilogx(frequencies, abs(squeeze(resp6))); 
grid on;
xlabel('Frequency');
ylabel('Magnitude');
title('Frequency Response');
legend("\xi_1 = 0.1","\xi_2 = 0.2","\xi_3 = 0.4","\xi_4 = 0.7","\xi_5 = 1", "\xi_6 = 1.5")
xlim([0,4])
hold off

figure_6.png