Too近期学习

qq_52501378

已于 2023-09-01 22:31:31 修改

阅读量315

点赞数 1

文章标签：学习

于 2023-08-29 21:15:30 首次发布

本文链接：https://blog.csdn.net/qq_52501378/article/details/132058920

版权

时域 ↔ 频域

傅里叶变换：

$F(\omega )=\int_{-\infty}^{+\infty}f(t)e_{}^{-j\omega t}dt$

逆傅里叶变换：

$f(t)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega )e^{j\omega t}d\omega$

我们可以对符合一定条件的连续函数进行傅里叶变换，但是很多时候我们获得的信号是离散的序列

离散傅里叶变换：

设 $x(n)$ 为有限长序列，只在0≤n<N-1时有值，现在把他扩充为周期为N的周期性序列 $\widetilde{x(n)}$

$X(k)=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}kn} (k=0,1,2,...,N-1)$

离散傅里叶反变换：

$x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k)e^{j\frac{2\pi }{N}nk} (0,1,2,...,N-1)$

matlab代码实现

利用快速傅立叶(FFT)的办法实现离散傅里叶(DFT)变换

Y = fft(X,n)

对输入信号X进行快速傅里叶变换

若未指定n，则返回值Y的长度将与X相同

指定n后，（将对X进行截断或补零到n位），返回值Y的长度与n相同

fs=8000;%采样频率
N=4000;%采样点数
t=0:1/fs:(N-1)/fs;%采样时间
f1=1000;f2=1010;%信号的频率为1000,1010hz
S=5*sin(2*pi*f2*t)+2+cos(2*pi*f1*t+pi/2);%给两个正弦信号
 
subplot(311)
plot(t,S)
 
%%fft
n=4000 %fft点数，一般应该是2的幂次
Y=fftshift(fft(S,n)); %把fft的零频分量移到屏幕中心
fshift=(-n/2:n/2-1)*fs/n;
%删除小幅值变换值,避免相位振荡
tol = 1e-6;
Y(abs(Y)<=tol)=0;

P1=abs(Y/n);
P2=2*P1;
P2(n/2+1)=P2(n/2+1)/2;
rad=angle(Y);
pha=rad/pi;

subplot(312)
plot(fshift,P2)
ylabel('Amplitude');
axis([950,1050,-inf,inf])
subplot(313)
stem(fshift,pha,'MarkerSize',1)
ylabel('Phase (\pi)');
axis([950,1050,-inf,inf])

采样频率： $f_s$

采样点数： $N$

fft的分辨率： $f_s/N$

实际运用时可以在原信号之后补零，从而改善分辨率(程度有限)

时域↔s域

电路分析中，我们经常需要求解线性微分方程，在实数域中进行这样的求解并不容易。为了计算方便，我们可以对实变量函数进行拉普拉斯变换，再在复数域中进行各种计算，最后所得的结果再逆变换到复数域中。

拉普拉斯变换：

$F(s)=L[f(t)]=\int_{0}^{+\infty}f(t)e^{-st}dt$

拉普拉斯逆变换：

$f(t)=L^{-1}[F(s)]=\frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}F(s)e^{st}ds$

其中σ是一个使F(s)的积分路径在收敛域内的实数

图源https://www.zhihu.com/search?q=%E6%8B%89%E6%99%AE%E6%8B%89%E6%96%AF%E9%80%86%E5%8F%98%E6%8D%A2&utm_content=search_history&type=content 拉普拉斯变换的微分性质：

$L[\frac{df(t)}{dt}]=\int_{0^-}^{+\infty}\frac{df(t)}{dt}e^{-st}dt=\int_{0^-}^{+\infty}e^{-st}df(t)=sF(s)-f(0^-)$

这里拉普拉斯积分的范围 $\int_{0^-}^{+\infty}$ ，从而确保脉冲函数 $\delta(t)$ 包含在积分内

$\int_{0}^{+\infty}e^{-st}df(t)=e^{-st}f(t)|_{0}^{+\infty}+s\int_{0}^{+\infty}f(t)e^{-st}dt$

在 $t\rightarrow +\infty$ 时， $e^{-st}f(t)=0$ ，故 $L[f^{'}(t)]=-f(0)+sL[f(t)]$

推广到n阶导， $L[f^{(n)}(t)]=-f^{(n-1)}(0)+sL[f^{n-1}(t)]$

可以进行多次迭代， $L[f^{(n)}(t)]=s^{n}L[f(t)]-f^{(n-1)}(0)-sf^{(n-2)}(0)-...-s^{n-1}f(0)$

如果 $f(t)$ 以及它的各阶导在t=0时都取零，则 $L[f^{(n)}(t)]=s^{n}L[f(t)]$

拉普拉斯变换的积分性质：

$L[\int_{0}^{t}f(t')dt']=\int_{0}^{+\infty}[\int_{0}^{t}f(t')dt']e^{-st}dt=-\frac{1}{s}\int_{0}^{+\infty}[\int_{0}^{t}f(t')dt']de^{-st}$

进行分部积分：

$\begin{aligned} L[\int_{0}^{t}f(t')dt']&=-\frac{1}{s}\int_{0}^{+\infty}[\int_{0}^{t}f(t')dt']de^{-st}\\&=-\frac{1}{s}\left \{ [(\int_{0}^{t}f(t')dt')e^{-st}]|_{0}^{+\infty} - \int_{0}^{+\infty}e^{-st}f(t)dt \right \}\\&=+\frac{1}{s}L[f(t)] \end{aligned}$

推广到n重积分：

$L[\underbrace{\int_{0}^{t}\int_{0}^{t}...\int_{0}^{t}}_{n \text{ times}}f(t)dt]=\frac{1}{s^{n}}L[f(t)]$

在电路分析中的应用：

图源：电子电路分析-拉普拉斯变换（应用篇） - 知乎

$\left\{\begin{matrix} E=&R_1i_1+L(\frac{di_1}{dt}-\frac{di_2}{dt}) \\\\0=&R_2i_2-L(\frac{di_1}{dt}-\frac{di_2}{dt}) \end{matrix}\right.$

对方程两边作拉普拉斯变换：

$\left\{\begin{matrix} \frac{E}{s}=R_1I_1(s)+L(sI_1(s)-i_1(0))-L(sI_2(s)-i_2(0)) \\ \\ 0=R_2I_2(s)-L(sI_1(s)-i_1(0))+L(sI_2(s)-i_2(0)) \end{matrix}\right.$

初始时刻电信号为0

$\left\{\begin{matrix} \begin{aligned}I_1(s)&=(\frac{E}{R_1+R_2})\frac{s+\frac{R_2}{L}}{s(s+\frac{R_1R_2}{L(R_1+R_2)})}\\\\I_2(s)&=(\frac{E}{R_1+R_2})\frac{1}{s+\frac{R_1R_2}{L(R_1+R_2)}}\end{aligned} \end{matrix}\right.$

进行拉普拉斯逆变换，得

$\left\{\begin{matrix}\begin{aligned} i_1(t)&=\frac{E}{R_1}-\frac{E}{R_1}(\frac{R_2}{R_1+R_2})e^{-\frac{R_1R_2}{L(R_1+R_2)}t} \\ \\ i_2(t)&=(\frac{E}{R_1+R_2})e^{-\frac{R_1R_2}{L(R_1+R_2)}t}\end{aligned} \end{matrix}\right.$

在初始条件不为零时求解线性微分方程：

来自：系统传递函数 - 搜索结果 - 知乎

线性微分方程 $\frac{d^2y}{dt^2}+5\frac{dy}{dt}+6y=6,\frac{dy}{dt}|_{t=0}=y(0)=2$

对微分方程两边进行拉普拉斯变换：

$s^2Y(s)-sy(0)-\frac{dy}{dt}|_{t=0}+5sY(s)-5y(0)+6Y(s)=\frac{6}{s}$

带入初始条件：

$Y(s)=\frac{2s^2+12s+6}{s(s^2+5s+6)}=\frac{2s^2+12s+6}{s(s+2)(s+3)}=\frac{1}{s}-\frac{4}{s+3}+\frac{5}{s+2}$

进行拉普拉斯反变换：

$y(t)=1-4e^{-3t}+5e^{-2t}(t>0)$

连续系统↔离散系统

原文：如何理解离散傅里叶变换及Z变换 - 知乎

模拟信号是指利用连续变化的物理量表示的信息，它在时域上是连续的；在计算机中，数字信号的大小常用有限位的二进制来表示，数字信号在时域上是离散的。

模拟信号：

数字信号：

离散时间傅里叶变换（DTFT）：

傅里叶变换是用来处理连续系统的，对于离散系统，我们需要进行离散时间傅里叶变换。假如 $x(t)$ 为连续系统信号，要转换成离散系统，需要先进行采样，采样频率为 $f_s$ ，冲击采样序列为：

$\delta_s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nT_s)$

则取样之后的信号为：

$x_s(t)=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s)$

连续信号的傅里叶变换定义为：

$X(j\omega )=\int_{-\infty}^{+\infty}x(t)e^{-j\omega t}dt$

则采样后的信号傅里叶变换为：

$X_s(j\omega)=\int_{-\infty}^{\infty}(\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT_s))e^{-j\omega t}dt$

交换一下积分与求和的顺序：

$X_s(j\omega)=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\delta(t-nT_s)(x(t)e^{-j\omega t})dt$

由δ(t)的筛选特性可知：

$X_s(j\omega)=\sum_{n=-\infty}^{\infty}x(nT_s)e^{-j\omega nT_s}$

对于数字信号而言，我们只有 $x(t)$ 采样后的信号 $\left \{ x[n] \right \}$ ，第n各采样发生在时间 $t=nT_s$ 。因此可以将连续信号的傅里叶变换写成如下形式：

$X_s(j\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-j\omega nT_s}$

上式就称为DTFT，也就是无限长离散信号的傅里叶变换。但是DTFT的结果频率ω还是连续的，不能直接应用在数字系统中。

离散傅里叶变换（DFT）

如果一个信号的频谱（DTFT）是连续的，我们仍然没法在机器中处理，这时就需要把频率也离散化，所需要的技术就是离散傅里叶变换DFT，即具有周期特性离散信号的傅里叶级数（就是将无限长的离散信号进行截断至N个采样点，然后将这N个采样点进行周期延拓），DFT与IDFT表达式如下：

$\begin{aligned} X[k]&=\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\frac{ 2\pi }{N}nk},&0\leqslant k< N-1\\ x[n]&=\sum_{n=0}^{N-1}X[k]e^{j\frac{2\pi}{N}nk},&0\leqslant n< N-1 \end{aligned}$

当然更多的时候也可以表示成：

$\begin{aligned}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}nk},&0\leqslant k< N-1\\ x[n]&=\frac{1}{N}\sum_{n=0}^{N-1}X[k]e^{j\frac{2\pi}{N}nk},&0\leqslant n< N-1 \end{aligned}$

即1/N可以放在DFT或者IDFT里面，两者是完全等价的。

这里从连续周期信号的傅里叶级数出发再次推导DFT，对于连续信号 $f(t)$ ，它的傅里叶级数：

$c_n=\frac{1}{T}\int_{0}^{T}f(t)e^{-j\frac{2\pi}{T}nt}dt$

对于连续信号 $x(t)$ 按照采样时间 $T_s$ 进行抽样N次，并将这N个数值进行周期延拓，可以得到周期的离散信号，在一个周期内，其表达式为：

$x_s(t)=x(t)\sum_{0}^{N-1}\delta(t-nT_s)$

可以获得离散周期信号的傅里叶级数为：

$X[jk\omega]=\frac{1}{T}\int_{-T/2}^{T/2}x(t)\sum_{n=0}^{N-1}\delta(t-nT_s)e^{-j\frac{2\pi}{T}kt}dt$

将积分与求和调换一下顺序：

$X[jk\omega]=\frac{1}{T}\sum_{n=0}^{N-1}\int_{-T/2}^{T/2}x(t)\delta(t-nT_s)e^{-j\frac{2\pi}{T}kt}dt$

δ函数具有筛选性质，上式可以离散化，我们有 $T\rightarrow NT_s,t\rightarrow nT_s$ ，上式可以简化为：

$X[jk\omega]=\frac{1}{NT_s}\sum_{n=0}^{N-1}x(nT_s)e^{-j\frac{2\pi}{NT_s}knT_s}=\frac{1}{NT_s}\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}$

令

$X[jk\omega]\cdot T_s=X[k]$

所以可以看到：

$X[k]=\frac{1}{N}\sum_{n=0}^{N-1}x(nT_s)e^{-j\frac{2\pi}{NT_s}knT_s}=\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}kn}$

这就是离散周期信号的傅里叶变换，理论上离散周期信号的频谱有无穷多，但是由于 $\left \{ e^{j\frac{2\pi}{N}kn} \right \}$ 的周期性，一般我们只取主值区间 $0\leqslant k< N-1$ 进行研究。

上面我们计算了离散周期信号的频谱，即不同频率下的系数，下面给出离散数值 $k$ 对应的频率：

$k\omega=k\frac{2\pi}{T}=k\frac{2\pi}{NT_s}=2\pi\frac{k}{N}\frac{1}{T_s}=2\pi f_s\frac{k}{N}$

也就是说：

$f_k=\frac{k}{N}f_s$

离散周期傅里叶变换（DFT）后的第k个数对应的频率是 $\frac{k}{N}f_s$ ，这个结论在后面Z变换的性质中会用到。快速傅里叶变换（FFT）和DFT本质上是一回事，只是它利用了DFT内部蕴藏的规律进行了一种简化的算法。

拉普拉斯变换和Z变换的关系

我们试着从收敛域的角度理解拉普拉斯变换和z变换的关系

信号 $f(t)$ 要保证能进行傅里叶变换，需要满足绝对可积条件：

$\int_{-\infty}^{\infty} \left | f(t) \right |< +\infty$

我们给 $f(t)$ 强制乘上一个衰减函数 $e^{-\sigma t}$ ，同时限制积分域从0开始，确保 $e^{-\sigma t}<1$ ，从而保证了：

$\int_0^{\infty} \left | f(t)e^{-\sigma t} \right |< +\infty$

此时 $\int_0^{\infty} f(t)e^{-\sigma t-j\omega t}=L[f(t)]$ ，

具体可见：从另一个角度看拉普拉斯变换 - 知乎

对于离散数列，我们试着给 $x[n]$ 乘一个指数函数 $e^{-\sigma\cdot nT_s}$ ，对应连续信号的衰减函数 $e^{-\sigma t}$ ，即用 $nT_s$ 对应时间 $t$ 。

则函数 $x[n]e^{-\sigma nT_s}$ 的DTFT为：

$X(j\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-\sigma nT_s}e^{-j\omega nT_s}$ ，

进一步化简：

$X(j\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-(\sigma+j\omega) nT_s}$

这看起来比较复杂，试着让它变得简洁，我们令

$z=e^{(\sigma+j\omega)T_s}$

则DTFT的表达式转变为双边Z变换的公式，我们可以发现 $\sigma=0$ 时，Z变换的表达式也就是离散时间傅里叶变换：

$X(j\omega )=\sum_{n=-\infty}^{\infty}x[n]z^{-n}$

我们知道在s变换中：

$s=\sigma+j\omega$

我们也知道：

$z=e^{(\sigma+j\omega)T_s}$

显然能得出以下结论：

$z=e^{sT_s}$

也就是说，z域是s域的进一步映射，z变换和拉普拉斯变换通过某种映射连接在一起。同时，由 $z=e^{sT_s}$ 还可以看出， $\sigma=0$ 时， $s=j\omega$ ，则 $z=e^{j\omega T_s}$ ，这时的z就是一个单位圆。

在深入研究一下z和s之间的映射关系：

由 $s=\sigma+j\omega$ ，当 $\sigma=0$ 时， $s=j\omega$ ，对应的是 $s$ 域的虚轴，而此时 $z=e^{j\omega T_s}$ 对应的是单位圆，也就是说z变换将s域的虚轴映射成z域的单位圆。

当 $\sigma >0$ 时， $s=\sigma+j\omega$ ，对应的是s域的正半轴，而此时 $z=e^{\sigma T_s}e^{j\omega T_s}$ ，由于 $e^{\sigma T_s}> 1$ ，也就是说此时z变换将s域的正半轴映射到了z域的单位圆外部。

当 $\sigma < 0$ 时， $s=\sigma+j\omega$ ，对应的是s域的负半轴，而此时 $z=e^{\sigma T_s}e^{j\omega T_s}$ ，由于 $e^{\sigma T_s}<1$ ，也就是说此时z变换将s域的负半轴映射到了z域的单位圆内部。

那么将s域映射到z域的好处是什么呢？

如果 $x(t)$ 和 $X(s)$ 是拉普拉斯变换对：

$X(s)=L[x(t)]$

则根据s变换的延时特性，可以得到：

$e^{-sT_s}X(s)=L[x(t-T_s)]$

现在 $x(t)$ 延时了 $n\cdot T_s$ ，那么：

$(e^{-sT_s})^{n}X(s)=L[x(t-n\cdot T_s)]$

可见，在时域进行延时操作，就相当于在频域做旋转操作， $e^{-sT_s}$ 大量出现在具有延时特性系统的拉普拉斯变换里面，如果用 $z=e^{sT_s}$ 来表示，计算会方便很多。简单来说，z就相当于定义了一个旋转操作符，这对于具有延时采样的离散系统特别有用。

系统函数

连续系统的系统函数

设系统的输入量为 $x_r(t)$ ，输出量为 $x_c(t)$ ，则它的传递函数是指初始条件为零时，输出量的拉普拉斯变换式 $X_c(s)$ 和输入量的拉普拉斯变换式 $X_r(s)$ 之比，并记作 $G(s)$ ：

$G(s)=\frac{L[x_c(t)]}{L[x_r(t)]}=\frac{X_c(s)}{X_r(s)}$

只要能写出系统的微分表达式，就可以求出其传递函数，例如常见的弹簧-质量-阻尼机械系统，令x(t)为位移，f(t)为质量块上的外力：

$M\frac{d^2x(t)}{dt^2}+R\frac{dx(t)}{dt}+Kx(t)=f(t)$

进行拉普拉斯变换：

$(Ms^2+Rs+K)X(s)=F(s)$

则系统的传递函数为：

$G(s)=\frac{1}{Ms^2+Rs+K}$

系统传递函数的极零点以及它的稳定性：

原文见：（4）传递函数零极点 - 知乎

一个有理传递函数，可表示为两个关于s的多项式之比：

$G(s)=\frac{b_1s^m+b_2s^{m-1}+...+b_{m+1}}{s^n+a_1s^{n-1}+...+a_n}$

可以通过因式分解转换成零极点的形式：

$G(s)=\frac{N(s)}{D(s)}=K\frac{\prod_{i=1}^{m}(s-z_i)}{\prod_{i=1}^{n}(s-p_i)}$

极点位于s平面中传递函数幅值无穷大处，即：

$\lim\limits_{s\to p_i}G(s)=\infty$

零点位于s平面中传递函数幅值为零处，即：

$\lim\limits_{s \to z_i}G(s)=0$

现在我们考虑单位脉冲信号δ(t)，有：

$L[\delta(t)]=1$

因此在单位脉冲信号输入下，系统输入Y（s），即系统的单位冲激响应：

$Y(s)=G(s)U(s)=G(s)$

因为单位冲激响应与系统传递函数具有相同的复域表达式，因此通过分析系统传递函数即可得到许多与系统动态表现直接相关观点结论。现在对于极点pi各不相同的清况，Y(s)可转化成分式相加的形式：

$Y(s)=\sum_{i=1}^{n}\frac{C_i}{s-p_i}$

为求解系数Ci，以C1为例：我们将上式左右同时乘s-p1，有：

$(s-p_1)Y(s)=C_1+\frac{(s-p_1)C_2}{s-p_2}+...+\frac{(s-p_1)C_n}{s-p_n}$

令s=p1，得：
$C_1=(s-p_1)Y(s)|_{s=p_1}$

以此类推，有：

$C_i=(s-p_i)Y(s)|_{s=p_i}$

确定系数后，再根据：

$L^{-1}[\frac{1}{s+a}]=e^{-at}$

便可以求得传递函数的拉普拉斯逆变换结果：

$\begin{aligned} y(t)&=L^{-1}[Y(s)]\\&=L^{-1}[\sum_{i=1}^{n}\frac{C_i}{s-p_i}]\\&=\sum_{i=1}^{n}C_ie^{p_it} \end{aligned}$

可以看到，指数部分仅仅由极点pi决定，即时间t的系数恰好是极点pi。因此，极点决定系统的动态表现与稳定性。

下面进行分类讨论，当极点为实数时，即位于s平面的实轴上。

①当 $p_i>0$ ，极点位于复平面右侧，有：

$\lim\limits_{t\to \infty}C_ie^{p_it}=\infty$

随着时间趋于无穷，发散速度由 $p_i$ 决定，当传递函数的所有极点中中出现一个这样的极点，则这个系统不稳定。

②当 $p_i<0$ ，极点位于复平面左侧，有：

$\lim\limits_{t\to \infty}C_ie^{p_it}=0$

其会随着时间收敛到0，收敛速度由 $p_i$ 决定，其绝对值越小（越靠近虚轴），收敛速度越慢。因此对于最终收敛的时域表达式来说，收敛的速度主要取决于那些更靠近0的极点。距离虚轴最近的极点所对应的响应分量在系统响应中起主导作用，因此称之为主导极点。

③当 $p_i=0$ ，极点位于原点，有：

$\lim\limits_{t\to \infty}C_ie^{p_it}=C_i$

表现为常值。

当极点的虚部不为零时，极点在复平面的实轴以外的位置。对于关于s的实系数多项式来说，其解一定为实数或共轭复数，因此有虚部的极点均为共轭极点，即关于实轴对称。设：

$\left\{\begin{matrix}\begin{aligned} p_1&=\sigma_0+j\omega_0 \\ p_2&=\sigma_0-j\omega_0\end{aligned} \end{matrix}\right.$

其单位冲激响应中关于该共轭极点的分量 $y_{1,2}(t)$ 为：

$\begin{aligned} y_{1,2}(t)&=C_1e^{(\sigma_0+j\omega_0)t}+C_2e^{(\sigma_0-j\omega_0)t})\\ &=e^{\sigma_0t}(C_1e^{j\omega_0t}+C_2e^{-j\omega_{0} t}) \\ &=e^{\sigma_0t}(C_1(cos(\omega_0t)+jsin(\omega_0t))+C_2(cos(\omega_0t)-jsin(\omega_0t)))\\ &=e^{\sigma_0t}((C_1+C_2)cos(\omega_0t)+(C_1-C_2)jsin(\omega_0t)) \end{aligned}$

计算C1和C2：(此处的公式我仍未看懂)

$\begin{aligned} C_1&=(s-p_1)Y(s)|_{s=p_1}\\ &=(s-p_1)Y^*(s)\frac{1}{(s-p_1)(s-p_2)}|_{s=p_1}\\ &=Y^*(s)\frac{1}{s-p_2}|_{s=p_1}\\ &=Y^*(p_1)\frac{1}{2j\omega_0} \end{aligned}$

$\begin{aligned} C_2&=(s-p_2)Y(s)|_{s=p_2}\\ &=(s-p_2)Y^*(s)\frac{1}{(s-p_1)(s-p_2)}|_{s=p_2}\\ &=Y^*(s)\frac{1}{s-p_1}|_{s=p_2}\\ &=Y^*(p_2)\frac{1}{-2j\omega_0} \end{aligned}$

其中Y*(p1)与Y*(p2)互为共轭，因此C1与C2互为共轭，有：

$C_1=C_r+jC_i,C_2=C_r-jC_i$

进一步的，有：

$\begin{aligned} y_{1,2}(t)&=e^{\sigma_0t}((C_1+C_2)cos(\omega_0 t)+(C_1-C_2)jsin(\omega_0t))\\ &=e^{\sigma_0t}(2C_rcos(\omega_0t)-2C_isin(\omega_0t))\\ &=C'e^{\sigma_0t}sin(\omega_0t+\phi ) \end{aligned}$

最终得到一个震荡的表达式，可以看出，极点的虚部意味着时域的振荡，同时决定振荡频率，实部收敛则决定收敛性。当 $\sigma_0<0$ 时分量振荡收敛，当 $\sigma_0>0$ 时响应分量振荡发散，当 $\sigma_0=0$ 时，响应分量稳定振荡。

极点分布情况对应的时域函数可总结为下图：

零点与输出响应：

考虑以下两个具有相同极点的系统 $G_1(s),G_2(s)$ ，此外他们还具有相同的稳态增益 $G_1(0)=G_2(0)=1$

$\begin{aligned} G_1(s)&=\frac{2}{(s+1)(s+2)}=\frac{2}{s+1}-\frac{2}{s+2}\\ G_2(s)&=\frac{2(s+1.1)}{1.1(s+1)(s+2)}=\frac{0.18}{s+1}+\frac{1.64}{s+2} \end{aligned}$

G2在极点 $p_1=-1$ 附近存在一零点 $z_1=-1.1$ 。通过 $s+1$ 分式项系数发现，该零点削弱了极点 $p_1=-1$ 在响应中的分量，一定程度起到抵消极点的作用，若该零点与极点相同则可完全抵消极点。

考虑更一般的情况，已知：

$C_i=(s-p_i)Y(s)|_{s=p_{i}}$

若存在零点位于 $p_i$ 附近，则 $Y(s)$ 的值会由于该零点变小，零点与 $p_i$ 越接近则影响越显著。零点除了在极点附近可以起抵消的作用，零点的加入也对系统阶跃响应有显著影响，考虑系统 $G(s)$ ：

$G(s)=\frac{24(s+z)}{z(s+4)(s+6)}$

可以将其拆成两个系统并联：

$G(s)=G_0(s)+G_1(s)=\frac{24}{(s+4)(s+6)}+\frac{24s}{z(s+4)(s+6)}$

其中 $G_0$ 与 $G$ 具有相同的极点与稳态增益，唯一的区别在是否存在零点，计算其阶跃响应：

动态性能：阶跃输入对系统来说是最严峻的工作状态。如果对于阶跃输入，系统的动态性能满足要求，那么系统在其他形式输入时，动态性能通常仍能满足要求

$\begin{aligned} Y(s)&=\frac{1}{s}G(s)\\ &=\frac{1}{s}(G_0(s)+G_1(s))\\ &=\frac{24}{s(s+4)(s+6)}+\frac{24}{z(s+4)(s+6)} \end{aligned}$

经过拉普拉斯逆变换可得：

$\begin{aligned} y_0(t)&=1-3e^{-4t}+2e^{-6t}\\ y_1(t)&=\frac{12}{z}e^{-4t}-\frac{12}{z}e^{-6t} \end{aligned}$

进一步的，有：

$y(t)=y_0(t)+y_1(t)=1+(\frac{12}{z}-3)e^{-4t}+(2-\frac{12}{z})e^{-6t}$

在z分别取2，3，4，5，6，36时，其阶跃函数的系统响应为：

可以看到，零点的加入可能使得本不包含共轭复极点的系统阶跃响应产生超调，随着极点距离原点的距离越近影响越显著，当零点位于复平面左侧无穷远时，对系统的影响无穷小。

奈奎斯特稳定判据

FB麦克风滤波器设计以及闭环系统的稳定性

图源：控制工程基础几个传递函数小总结以及一些汇总 - 知乎

参考：奈奎斯特稳定性判据简单清晰的理解和推导 - 知乎

在控制系统理论中，开环传递函数(open loop gain)指的是： $G(s)H(s)$

包含负反馈的闭环传递函数： $\frac{G(s)}{1+H(s)G(s)}$

线性时不变系统稳定的充分必要条件是：传递函数的极点位于负半平面（闭环传递函数极点的实部小于零）。这也是奈奎斯特稳定性判据的理论基础。

传递函数的极点都在处在负平面时，系统是稳定的；现在我们想知道闭环传递函数的稳定性，直观上只需要观察 $\frac{G(s)}{1+H(s)G(s)}$ 的极点位置就可以。但是也可以从开环传递函数 $G(s)H(s)$ 出发，利用奈奎斯特稳定判据判断系统的稳定性。具体实施见参考文。

根据理论分析，我们已经得到N=Z-P这个关系。N值可以通过奈奎斯特曲线绕（-1，j0）点的情况获得。P值也可以通过观察开环传函表达式中极点实部大于零的个数获得。这样就可以得到Z值，Z=N+P。而Z值就是闭环传函极点的实部大于零的个数。当Z得0时，则闭环传函极点实部都小于零，系统稳定。若Z大于零，闭环传函有极点实部大于零，系统不稳定。这正是奈奎斯特稳定判据的含义

裕度判据

参考：

增益裕度和相位裕度的理解 - 知乎

带你搞懂自控频域中的稳定性分析 - 知乎

“系统设计的相位裕度太低了，容易振荡”，“相位裕度太低了，参数振动时，系统容易不稳定”。我们试着从闭环系统出发理解这句话。

闭环系统的结构，X(s)为输入，Y(s)为输出。

$[X(s)-Y(s)*H(s)]*G(s)=Y(s)$

$\frac{Y(s)}{X(s)}=\frac{G(s)}{1+H(s)G(s)}$ ，这里为了分析耳机内的FB降噪，令G(s)=1

首先我们需要理解相位对负反馈的影响，我们令H(s)附加的相位延迟为-180°，此时我们给定一个正弦信号输入，信号经过H(s）后再次反相叠加到原信号上。这时可以发现我们的系统实际上成为了一个正反馈系统，其中H(s)的绝对值就决定了正反馈的大小;正反馈的系统不一定是发散的，当H(s)>1时，系统不稳定；H(s)=1时，系统临界稳定。

一般我们会用dB来形容H(s)的增益，临界稳定时， $20lg\left |H(s)\right |=0dB$

相位裕度定义：

$20lg\left |H(s)\right |=0dB$ 时，对应的角频率为 $\omega_c$ ，此时的相位响应距离-180°的角度， $\displaystyle PM=\phi (\omega_c )-(-180^{\circ})=180^{\circ}+\phi (\omega_c )$

增益(幅值)裕度定义：

相位响应为-180°时，此时的角频率对应的幅度响应距离0dB的大小

$GM=-[20lg\left|H(s)\right|]$

耳机的反馈控制器

耳机FB控制器优化，利用负反馈控制进行降噪，此时拾取的是FB处的噪声。降噪的目的是消除耳膜(Ear)处的噪声，但实际上我们的负反馈通路拾取到的信号点FB距离耳膜还有一段距离，因此如何补偿两点的差异成了降噪工作的重点。

$Noise(s)*Noise2FB(s)-Y_1(s)*Filter(s)*Spk2FB(s)=Y_1(s)$

$\frac{Y_1(s)}{Noise(s)}=\frac{Noise2FB(s)}{[1+Filter(s)*Spk2FB(s)]}$

为了与上面的简化系统类比，我们令 $\displaystyle X(s)=Noise(s)*Noise2FB(s)$

$G(s)=1,H(s)=Filter(s)*Spk2FB(s)$

因此我们可以根据 $H(s)=Filter(s)*Spk2FB(s)$ 来判断FB系统是否会导致啸叫

耳机的反馈控制器优化

由于FB点的位置和人耳鼓膜点不同，我们需要尽可能缩小二者差异

对于FB处的信号Y1，我们有

$\displaystyle Noise(s)*Noise2FB(s)-Y_1(s)*Filter(s)*Spk2FB(s)=Y_1(s)$

Y1与经过被动降噪后Noise的比值

$\displaystyle \frac{Y_1(s)}{Noise(s)*Noise2FB(s)}=\frac{1}{1+Filter(s)*Spk2FB(s)}$

对于Ear处的信号Y2，由于无法采集到Ear点的噪声，在反馈回路中我们采用Y1替代Y2

$\displaystyle Noise(s)*Noise2Ear(s)-Y_1(s)*Filter(s)*Spk2Ear(s)=Y_2(s)$

Y2与经过被动降噪后Noise的比值

$\displaystyle \frac{Y_2(s)}{Noise(s)*Noise2Ear(s)}=1-\frac{Y_1(s)*Filter(s)*Spk2Ear(s)}{Noise(s)*Noise2Ear(s)}$

将Y1用FB回路中的等式替代

$\begin{aligned}\frac{Y_2(s)}{Noise(s)*Noise2Ear(s)}&=1-\frac{Filter(s)*Spk2Ear(s)}{Noise(s)*Noise2Ear(s)}*\frac{Noise(s)*Noise2FB(s)}{1+Filter(s)*Spk2FB(s)}\\&=1-\frac{Filter(s)*Spk2Ear(s)}{Noise2Ear(s)}*\frac{Noise2FB(s)}{1+Filter(s)*Spk2FB(s)}\\&=\frac{1+Filter(s)*Spk2FB(s)-\frac{Filter(s)*Spk2Ear(s)*Noise2FB(s)}{Noise2Ear(s)}}{1+Filter(s)*Spk2FB(s)} \end{aligned}$

现在我们令 $\Delta(z)=\frac{Spk2Ear(s)*Noise2FB(s)}{Spk2FB(s)*Noise2Ear(s)}$ ，并且在 $Filter(s)*Spk2FB(s)>>1$ 的情况下

$\displaystyle \begin{aligned}\frac{Y_2(s)}{Noise(s)*Noise2Ear(s)}&=\frac{1+Filter(s)*Spk2FB(s)(1-\Delta(z))}{1+Filter(s)*Spk2FB(s)}\\&\approx \frac{Y_1(s)}{Noise(s)*Noise2FB(s)} +(1-\Delta(z))\end{aligned}$

此时可以看到 $\Delta(z)$ 趋向于1时，

$\displaystyle \frac{Y_2(s)}{Noise(s)*Noise2Ear(s)}\approx \frac{Y_1(s)}{Noise(s)*Noise2FB(s)}$

耳机的前馈控制器

耳机的前馈降噪，此时拾取的是外界的噪声。

其中Noise2FB为耳机的被动降噪PNC；FF-FR为前馈麦克风的频响；FF-filter是FF滤波器；Spk2FB是扬声器到FB的频响；开启前馈后耳内FB点的噪声为Y1。

$\begin{aligned} Y_1(s)&=Noise(s)*Noise2FB\\&+Noise(s)*(FF-FR)*(FF-Filter)*Spk2FB \end{aligned}$

我们希望耳内噪声尽可能小，因此可以得到FF目标滤波器，一般近似认为FF的频响是平直的

$\displaystyle FF-Filter=-\frac{Noise2FB}{(FF-FR)*(Spk2FB)}\approx -\frac{Noise2FB}{Spk2FB}$

根轨迹法

参考：

自动控制原理（四）_闭环特征方程_童年吹梦的博客-CSDN博客

【自控原理】第四章根轨迹法_二进制人工智能的博客-CSDN博客

根轨迹定义：根轨迹是指当系统开环的某个参数（如开环增益）从零变化到无穷大时，闭环特征方程的根在平面上移动的轨迹

离散系统的系统函数

参考：

第二章 z变换六、离散系统的系统函数、系统的频率响应_z变换求频率响应_Clark-dj的博客-CSDN博客

离散系统的系统函数及系统的频率响应（2-8）_序列频率响应_不尽木的博客-CSDN博客

信号与系统(Python) 学习笔记 (8.1) 离散系统z域分析 -- 系统函数 H(z)_Varalpha的博客-CSDN博客

在离散系统中，我们可以给出时域和z域的转换， $y[n]=x[n]*h[n],Y(z)=X(z)H(z)$

系统的单位脉冲响应h[n]的z变换可以写作：

$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\frac{Y(z)}{X(z)}$

现在考虑系统的频率响应 $H(e^{j\omega})$ ，即单位圆上的系统函数，或者h[n]的DTFT

$H(e^{j\omega})=H(z)|_{z=e^{j\omega}}=DTFT[h(n)]$

实际工作中，我们的信号往往是离散的，因此讨论离散系统的系统函数是很有必要的。

从连续系统的常系数微分方程出发：

$\sum_{k=0}^{N}a_k\frac{d^k y(t)}{dt^k}=\sum_{k=0}^Mb_k\frac{d^kx(t)}{dt^k}$

在离散系统中，我们给出它所对应的N阶(阶次是输出y的最高阶导数)线性常系数差分方程：

$\sum_{k=0}^Na_ky[n-k]=\sum_{m=0}^Mb_mx[n-m]$

取z变换得：

$\sum_{k=0}^Na_kz^{-k}Y(z)=\sum_{m=0}^Mb_mz^{-m}X(z)$

所以写出系统函数：

$\begin{aligned}H(z)&=\frac{Y(z)}{X(z)}\\&=\frac{\sum_{m=0}^Mb_mz^{-m}}{\sum_{k=0}^Na_kz^{-k}}\\&=K\frac{\prod_{m=1}^{M}(1-c_mz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})} \end{aligned}$

H(z)的零点：

$z=c_m,m=1,2,...,M$

H(z)的极点：

$z=d_k,k=1,2,...,N$

离散系统的H(z)的零极点与h(k)的关系：

极点在单位圆内（实轴上）：

$H(z)=\frac{Az}{z-a}\Rightarrow h[k]=Aa^ku[k],\left | a \right |<1$

极点在单位圆内（非实轴上），此时极点必共轭：

$H(z)=\frac{Az}{z-re^{j\beta}}+\frac{A^*z}{z-re^{-j\beta}}\Rightarrow 2\left | A \right |r^kcos(\beta k+\theta)u[k],r<1$

可以看到函数h[k]指数衰减， $\sum_{k=-\infty}^{+\infty} \left |h[k] \right |<\infty$ ，此时对于一个有界的输入x[n]，系统是稳定的。

如果有一个离散时间且有因果性的滤波器，其传递函数的极点落在复数z平面的单位圆内，此滤波器则为稳定的

采样定理

原文：如何理解离散傅里叶变换及Z变换 - 知乎

我们知道，两个信号在时域相乘，在频域相当于卷积；在时域卷积，在频域相当于相乘。

离散信号，其实就是连续信号 $x(t)$ 与狄拉克梳状函数（也就是采样函数）的相乘，这就是采样。采样是一个时域的相乘行为，那么在频域它就会对应卷积。

一个信号的频谱函数和狄拉克梳状函数卷积，结果会是什么？

——频谱函数被周期延拓了，而且延拓周期就是采样频率

举个例子：

对于余弦函数而言，比如 $x(t)=Acos(2\pi\cdot8\cdot t)$ ，其傅里叶变换包含两个频率分量，分别是8Hz和-8Hz，如下图：

我们用一个20Hz的狄拉克梳状函数(采样函数)对余弦信号进行抽样，采样函数为：

$\delta_s(t)=\sum_{n=0}^{N}\delta(t-nT_s)$ ，其中 $T_s=1/f_s ,f_s=20Hz$

其时域信号以及频谱为：

可见，抽样信号的频谱也呈现冲击信号的样貌，频谱中两个冲击串的间隔就是采样频率。

采样后的信号为：

$x_s(t)=x(t)\delta_s(t)=Acos(2\pi\cdot8\cdot t)\sum_{n=0}^{N}\delta(t-nT_s)$ ：

采样后的时域信号及频谱为：

可见，采样后的信号被周期延拓了，延拓的周期就是20Hz，也就是采样频率。

把三张图放在一起

通过以上分析，我们可以得出结论：对一个连续信号的采样，采样后的频谱相当于将采样前的频谱进行了延拓，延拓的周期就是采样频率

假设一个信号的频谱如下：

频谱中最大的频率为 $f_{max}$ ，用一个周期为 $f_s$ 的狄拉克梳状函数进行采样后的频谱为原频谱的周期延拓，示意图如下：

看起来很完美，但如果出现了一种情况 $f_{max}> \frac{1}{2}f_s$ ：

也就是说，原始频谱信号经过周期延拓后会有一部分重叠，这样在重叠的部分就会有信息的丢失，也就无法进行信号复原了，也就是说，对于连续信号的进行抽样离散的话，必须保证采样频率是原连续信号最大频率分量的2倍频率以上，否则信号就难以复原。这就是采样定理，又叫奈奎斯特采样定理或香农采样定理。

s域↔z域

原文：传递函数的离散化（以一阶低通RC滤波器为例）_传递函数离散化_nova_peng的博客-CSDN博客

在学习过程中，接触到的都是连续系统，而计算机是离散的，因此对于连续系统进行离散处理是很有必要的。描述连续系统使用的是微分方程，而描述离散系统采用差分方程，因此将系统离散化的一般流程为：

z变换是对离散序列的一种数学变换，常用来求解线性时不变的差分方程，它在离散系统中的地位等同于拉普拉斯变换在连续系统中的地位。

双边z变换：

$X(z)=Z\left \{ x[n] \right \}=\sum_{n=-\infty}^{+\infty}x[n]z^{-n}$

单边z变换：

$X(z)=Z\left \{ x[n] \right \}=\sum_{n=0}^{+\infty}x[n]z^{-n}$

逆z变换：

$x[n]=\frac{1}{2\pi j}\oint X(z)z^{n-1}dz$

原文：连续系统的离散化方法 - 知乎

1.反向差分法

对于给定面积计算，采用积分形式：

$u(t)=\int_{0_-}^{+\infty}e(\tau)d\tau$

我们给出：

$D(s)=\frac{U(s)}{E(s)}=\frac{1}{s}$

其微分方程为 $\frac{du(t)}{dt}=e(t)$ ，用反向差分替代微分，并对数据做离散化，可得：

$\frac{du(t)}{dt}\approx \frac{u[k]-u[k-1]}{T}=e[k]$

对上式两边取z变换得 $(1-\frac{1}{z})U(z)=TE(z)$ ，即

$D(z)=\frac{U(z)}{E(z)}=\frac{1}{\frac{1-z^{-1}}{T}}$

比较两式可得：

$s=\frac{1-z^{-1}}{T},D(z)=D(s) |_{s=\frac{1-z^{-1}}{T}}$

另外还可以将 $z^{-1}$ 作级数展开：

$z^{-1}=e^{-Ts}=1-Ts+\frac{T^2s^2}{2}-...$

取一阶近似 $z^{-1}\approx 1-Ts$ ，

S平面的稳定域可以通过 $s=\frac{1-z^{-1}}{T}$ 映射到z平面。因为s平面的稳定域为 $Re(s)<0$ ，因此可以看出z平面中的稳定域为：

$Re(\frac{1-z^{-1}}{T})=Re(\frac{z-1}{Tz})<0$

T为正数，将z写成 $z=\sigma+j\omega$ ，注意此时的 $\sigma$ 和 $\omega$ 不是s平面的实部和虚部，上式可以写成

$Re(\frac{\sigma+j\omega-1}{\sigma+j\omega})<0$

即

$Re[\frac{(\sigma+j\omega-1)(\sigma-j\omega)}{(\sigma+j\omega)(\sigma-j\omega)}]=Re[\frac{\sigma^2-\sigma+\omega^2+j\omega}{\sigma^2+\omega^2}]=\frac{\sigma^2-\sigma+\omega^2}{\sigma^2+\omega^2}<0$

上式可以写成

$(\sigma-\frac{1}{2})^2+\omega^2<(\frac{1}{2})$

由上式可以看出，s平面的稳定域映射到z平面上以 $\sigma=\frac{1}{2},\omega=0$ 为圆心，1/2为半径的圆内。

反向差分变换方法的主要特点如下：

1.变换计算简单

2.图1.1可以看出，s平面的左半平面映射到z平面的单位圆内部一个小圆内，因而，如果D(s)稳定，则变换后的D(z)也是稳定的

3.离散滤波器的过程特性及频率特性同原连续滤波器比较有一定的失真，需要较小的采样周期T

2.正向差分变换法

对于给定的

$D(s)=\frac{U(s)}{E(s)}=\frac{1}{s}$

其微分方程为 $\frac{du(t)}{dt}=e(t)$ ，用正向差分替代微分，并对数据做离散化，可得：

$\frac{du(t)}{dt}\approx \frac{u[k+1]-u[k]}{T}=e[k]$

两边取z变换得： $(z-1)U(z)=TE(z)$ ，即

$D(z)=\frac{U(z)}{E(z)}=\frac{1}{\frac{z-1}{T}}$

可以发现，在对 $D(s)$ 进行正向差分变换时，将其中的s直接用

$s=\frac{z-1}{T},D(z)=D(s)|_{s=\frac{z-1}{T}}$

另外还可以将z级数展开：

$z=e^{sT}=1+Ts+\frac{T^2s^2}{2}+...$

取一阶近似 $z\approx 1+Ts$ ，也可以得到 $s=\frac{z-1}{T}$ ，使用正向差分方法时，有个严重问题是，s平面的左半平面映射到z平面的单位圆外。因为s平面的稳定域为 $Re(s)<0$ ，这时z平面的稳定域为：

$Re(\frac{z-1}{T})<0$

令 $z=\sigma+j\omega$ ，则上式可以写成：

$Re(\frac{\sigma+j\omega-1}{T})<0$

也就是说， $\sigma<1$ ，此时的收敛域如图所示：

由此，得出正向差分法变换的特点：s平面左半平面的极点可能映射到z平面单位圆外。因而用这种方法所得到的离散滤波器可能是不稳定的，实际应用中基本上不采用这种方法。

3.双线性变换法

参考：

双线性变换（Tustin transform/bilinear transformation）_双线性变换法_teengad的博客-CSDN博客

双线性变换法详解 - 知乎

双线性变换（Tustin transform/bilinear transformation）_双线性变换法_teengad的博客-CSDN博客

双线性变换法，是一种基于梯形积分规则的数字积分变换方式，它能够防止频率混叠。

由梯形的面积近似积分可以得：

$u(kT)=u[(k-1)T]+\frac{T}{2}[e[(k-1)T]+e(kT)]$

其中 $u(kT)$ 为到kT时刻的阴影总面积， $e(kT)$ 为曲线对应的纵坐标值，T就是采样间距。对上式进行Z变换:

$U(z)-z^{-1}U(z)=\frac{T}{2}[E(z)+z^{-1}E(z)]$

双线性变换的本质就是用一阶估计法将连续时间传递函数D(s)中的s替换为离散值z，如下所示：

$D(z)=\frac{U(z)}{E(z)}=\frac{T}{2}\frac{1+z^{-1}}{1-z^{-1}}=D(s)|_{s=\frac{2(1-z^{-1})}{T(1+z^{-1})}}$

由此得到双线性变换公式：

$s=\frac{2(1-z^{-1})}{T(1+z^{-1})}$

另外还可以由Z变换定义 $z=e^{sT}$ 形式出发：

$e^{sT}=\frac{e^{\frac{Ts}{2}}}{e^{-\frac{Ts}{2}}}$

然后将分子和分母同时展开成泰勒级数，取前两项，也可以得到双线性变换公式（分子和分母都是s的线性表达）：

$z=\frac{1+{\frac{Ts}{2}}}{1-{\frac{Ts}{2}}}, \Rightarrow s=\frac{2(1-z^{-1})}{T(1+z^{-1})}$

在传递函数D(s)=1/s的稳定域中，传递函数极点小于零，可以得到如下关系：

$Re(\frac{2(1-z^{-1})}{T(1+z^{-1})})=Re(\frac{2}{T}\frac{z-1}{z+1})<0$

因为T＞0，上面不等式可以简化：

$Re(\frac{2}{T}\frac{z-1}{z+1})=Re(\frac{\sigma+j\omega-1}{\sigma+j\omega+1})=Re(\frac{\sigma^2-1+\omega^2+2j\omega}{(\sigma+1)^2+\omega^2})<0$

即 $\sigma^2+\omega^2<1$ ，稳定域图示如下：

可以看到，这相当于将s平面的负半区域映射到z平面的单位圆内部。

如果有一个连续时间且有因果性的滤波器，其传递函数D(s)的极点落在复数s平面的左半边，此滤波器则为稳定的。

如果有一个离散时间且有因果性的滤波器，其传递函数D(z)的极点落在复数z平面的单位圆内，此滤波器则为稳定的。

双线性变换将复数s平面的左半边映射到复数z平面的单位圆内，因此稳定的连续时间滤波器被转变成离散时间滤波器后也保有稳定性。

现在考虑离散系统的频率响应

$Y(z)=H(z)X(z)\Rightarrow Y(e^{j\omega })=H(e^{j\omega})X(e^{j\omega })$

令 $s=j\Omega$ ，T是采样周期，经过双线性变换：

$z=\frac{1-(\Omega T/2)^2}{1+(j\Omega T/2)^2}+j\frac{\Omega T}{1+(j\Omega T/2)^2}$

z的模为1，即s平面的虚数轴映射为z平面的单位圆。

设 $z=e^{j\omega},\omega$ 为归一化的数字圆频率， $\omega=2\pi \frac{f_d}{f_s}=2\pi {f_d}{T};f_d$ 为其归一化数字频率， $f_d=\frac{\omega}{2\pi T}$ ,T为采样间隔。

将z代入双线性变换公式：

$s=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}=j\frac{2}{T} tan\frac{\omega}{2}=\sigma+j\Omega=j\Omega$

于是：

$\Omega=\frac{2}{T}tan\frac{\omega}{2}$

也可以写为：

$\omega=2arctan(\frac{\Omega T}{2})$

上式中 $\omega$ 是数字角频率， $\Omega$ 是模拟角频率。

可以看出所有的数字角频率都被压缩到± $\pi$ 之内，但是却依然能保持单值对应,虽然消灭了频谱混叠，但是由于这个出发点是在原点进行泰勒展开做近似的，最后的频谱会存在一定的失真。

使用双线性变换将s平面的虚轴映射到z平面单位圆上时，模拟角频率变化到数字角频率

$\Omega\in (-\infty,\infty)\Rightarrow \omega\in(-\pi,\pi)$ 是非线性。这会造成频响峰值的频率发生畸变，所以要量化二者关系，用于畸变补偿。

补偿的基本思想：在D(s)变换成D(z)之前，将D(s)的断点频率预先加以修正（预畸变）使得修正后的D(s)变换成D(z)时正好达到所要求的断点频率。预畸变双线性变换的特点为：

以二阶系统为例，传递函数为H(s)，信号与系统第九章9.4.2

$m\frac{d^2y(t)}{dt^2}+R\frac{dy(t)}{dt}+ky(t)=x(t)\Rightarrow H(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_n s+\omega_n^2 }$

带入双线性变换 $s=\frac{2}{T}\frac{z-1}{z+1}=2f_s\frac{z-1}{z+1}$ ：

$H(z)=\frac{\omega_n^2}{[\frac{2}{T}\frac{z-1}{z+1}]^2+2\zeta\omega_n[\frac{z-1}{z+1}]+\omega_n^2}=\frac{\omega_n^2(z^2+2z+1)}{(4f_s^2+4\zeta\omega_n f_s+\omega_n^2)z^2+(-8f_s^2+2\omega_n^2)z+(4f_s^2-4\zeta \omega_n f_s+\omega_n^2)}$

对于畸变补偿，由于 $\Omega$ 和 $\omega$ 是耦合的，不能单独更改，所以引入系数K，使：

$\Omega_1=K\frac{2}{T}tan(\frac{\omega_1}{2})\leftrightharpoons K=\frac{T}{2}\frac{\Omega_1}{tan(\omega_1T/2)}$

那么畸变补偿的双线性变换为：

$s=\frac{\Omega_1}{tan(\frac{\omega_1T}{2})}\frac{z-1}{z+1}$

带入传递函数 $H(s)$ ，并令 $2M=\frac{\Omega_1}{tan(\omega_1T/2)}$ ：

$H(z)=\frac{\omega_n^2(z^2+2z+1)}{(4M^2+4\zeta\omega_n M+\omega_n^2)z^2+(-8M^2+2\omega_n^2)z+(4M^2-4\zeta \omega_n M+\omega_n^2)}$

设置参数 $\omega_n=2\pi*30; \zeta=0.1;\omega_1=2\pi*500$ ，画个图来看看差异

4.脉冲响应不变法

原文

https://isip.bit.edu.cn/docs/2016-05/20160511090607824812.pdf

基本原理：从时域响应出发，使离散系统的单位脉冲响应h[n]模仿连续系统的单位脉冲响应h(t)，h[n]等于h(t)的取样值。

（1）将连续系统函数H(s)进行部分分式展开

$H(s)=\sum_{k=1}^{N}\frac{A_k}{s-p_k}$

（2）对H(s)进行拉普拉斯反变换，得到系统的单位脉冲响应

$h(t)=\sum_{k=1}^NA_ke^{p_kt}u(t)$

（3）以时间间隔T对h(t)时域采样得到h(n)，即维持脉冲响应不变

$h[n]=\sum_{k=1}^NA_ke^{p_knT}u(nT)=\sum_{k=1}^NA_k(e^{p_kT})^nu(n)$

（4）对h(n)进行z变换

$\displaystyle H(z)=\sum_{k=1}^N\frac{A_k}{1-e^{p_kT}z^{-1}}$

得到变换关系：

$\displaystyle \frac{1}{s-s_k}\Leftrightarrow \frac{1}{1-e^{s_kT}z^{-1}}=\frac{z}{z-e^{s_kT\displaystyle }}$

可以看到s平面的极点映射到z平面

$p_k\rightarrow e^{p_kT}$

可以看到，有如下s平面到z平面的变换

$z=e^{sT}\Rightarrow e^{j\omega}=e^{j\Omega T}$ ，

模拟和数字频率由下式联系： $\omega=\Omega T$ 或 $e^{j\omega}=e^{j\Omega T}$

因此，一个数字频率可能有多个模拟频率等价，频谱有混叠

5.RC低通滤波器的离散形式：

原文：

传递函数的离散化（以一阶低通RC滤波器为例）_传递函数离散化_nova_peng的博客-CSDN博客

这里以一阶RC低通滤波器为例，得到连续系统到差分方程的转化。

首先给出低通滤波器的结构：

系统满足的微分方程： $x(t)=RC\frac{dy(t)}{dt}+y(t)$

进行拉普拉斯变换： $X(s)=(RCs+1)Y(s)$

传递函数： $H(s)=\frac{1}{1+RCs}=\frac{\frac{1}{RC}}{s+\frac{1}{RC}}$

记 $\frac{1}{RC}=a,H(s)=\frac{a}{s+a}$

离散化方法有很多，比如一阶向后差法、双线性变换法、零极点匹配法、保持器等价法、脉冲响应不变法、阶跃响应不变法。各个方法的特点的与优势不尽相同，下面将使用脉冲响应不变法、向后差分、双线性变换法进行系统的离散化。

使用反向差分法进行离散

连续系统的微分用采样时间为T的向后差分代替，并给出s域到z域的变换算子：

$\frac{dy(t)}{dt}=\frac{y[n]-y[n-1]}{T},s=\frac{1-z^{-1}}{T}$

将算子代入H(s),变换结果如下：

$H(z)=\frac{Y[z]}{X[z]}=\frac{a}{\frac{1-z^{-1}}{T}+a}=\frac{aT}{1+aT-z^{-1}}$

进行z反变换得差分方程：

$y[n]=\frac{aT}{1+aT}x[n]+\frac{1}{1+aT}y[n-1]$

双线性变换法

双线性变换法算子：

$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$

带入H(s）后的变换结果：

$H(z)=\frac{Y[z]}{X[z]}=\frac{a}{\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}+a}=\frac{\frac{aT}{aT+2}z+\frac{aT}{aT+2}}{z+\frac{aT-2}{aT+2}}$

进行反变换，解得：

$y[n]=\frac{aT}{aT+2}x[n]+\frac{aT}{aT+2}x[n-1]-\frac{aT-2}{aT+2}y[n-1]$

脉冲响应不变法

根据脉冲响应不变法的公式：

$H(z)=Z(H(s))$

$H(s)=\sum_{i=1}^{N}\frac{A_i}{S+S_i}\Rightarrow H(z)=\sum_{i=1}^{N}\frac{A_i}{1-e^{-s_iTz^{-1}}}$

可以得到系统传递函数的z变换如下：

$H(s)=\frac{a}{s+a},H(z)=\frac{Y[z]}{X[z]}=\frac{a}{1-e^{-aT}z^{-1}}$

进而得到差分方程：

$y[n]=ax[n]+e^{-aT}y[n-1]$

四种系统函数的matlab实现

画出连续传递函数和各种离散形式的效果

参考：

Matlab中lsim函数使用_lsim函数matlab_jk_101的博客-CSDN博客

代码：

clc;clear;close all;

figure;
a=0.3;T=1;
sys=tf(a,[1,a]);
sys
t=0:0.01:90;
xt=1*cos(0.3*t)+1*cos(2*t);
%lism(sys,u,t)绘制动态系统模型对输入历史记录(t,u)的模拟时间响应。
lsim(sys,xt,t);
hold on

n=0:(90/T);
xn=1*cos(0.3*n*T);
%反向差分法
%差分方程中x(n)项的系数
B=(a*T)/(1+a*T);
%差分方程中y(n)项的系数
A=[1,-1/(1+a*T)];
yn=filter(B,A,xn);
plot(n*T,yn,'Color','g','Marker','o');
%脉冲响应不变法
B=a;
A=[1,-exp(-a*T)];
yn2=filter(B,A,xn);
plot(n*T,yn2,'b+');
%双线性不变法
B=[(a*T)/(2+a*T),(a*T)/(2+a*T)];
A=[1,(a*T-2)/(a*T+2)];
yn3=filter(B,A,xn);
plot(n*T,yn3,'r*');


%标注
legend('tout','(y[n]-y[n-1])/T','imp','tustin')
%legend('(y[n]-y[n-1])/T','imp','tustin')
% 打开网格
grid on;

软磁材料

滤波器设计

阶跃函数在理想情况下可以被用作滤波器，因为它可以完美地分割频率。然而，在实际应用中，直接使用阶跃函数作为滤波器会遇到一些问题：

1. 非因果性：理想的阶跃函数滤波器是非因果的，也就是说，它需要未来的输入信息来计算当前的输出。在实际的信号处理中，我们无法获取未来的信息，因此无法实现非因果滤波器。

2. 无限冲激响应：理想的阶跃函数滤波器具有无限的冲激响应。这意味着，为了精确地实现阶跃函数滤波器，我们需要无限长的时间和无限多的计算资源。

3. 吉布斯现象：当试图用有限的频率成分来近似阶跃函数时，会出现所谓的吉布斯现象，即在阶跃点附近出现振荡和超调。

因此，虽然阶跃函数在理论上是一个理想的滤波器，但在实际应用中，我们通常使用其他类型的滤波器（如Butterworth滤波器，Chebyshev滤波器，椭圆滤波器等），这些滤波器在有限的资源下可以提供更好的性能。

原文：

数字信号处理|Matlab设计巴特沃斯低通滤波器（冲激响应不变法和双线性变换法）_双线性变换法设计低通滤波器_zzztutu的博客-CSDN博客

用双线性变换法设计IIR 滤波器MATLAB实现_timerring的博客-CSDN博客

https://isip.bit.edu.cn/docs/2016-05/20160511090607824812.pdf

设计流程
数字滤波器指标→模拟滤波器指标→合适的模拟滤波器→离散化成为数字滤波器

最小相位系统

参考：

数字滤波器(二)--最小相位延时系统和全通系统_全通最小相位系统零极点_爱吃骨头的猫、的博客-CSDN博客

当输入信号经过一个滤波器系统时，其相位可能会发生一定的移动，那么一个系统的相位到底跟什么有关，现在给定一个系统传递函数：

$\displaystyle H(z)=\frac{\sum_{m=0}^{M}b_rz^{-m}}{\sum_{k=0}^{N}a_kz^{-k}}=K\cdot z^{N-M}\frac{\prod_{m=1}^{M}(z-c_m)}{\prod_{k=1}^{N}(z-d_k )}$

可以看到它的相位响应：

$arg[\frac{H(e^{j\omega})}{K}]=\omega(N-M)+\sum_{m=1}^{M}arg(e^{j\omega}-c_m)-\sum_{k=1}^{N}arg(e^{j\omega}-d_k)$

可以看到H(z)共有N个极点和M个零点，所谓最小相位系统就是对输入信号的延时最小的系统。首先做一些定义。

$m_i$ ：单位圆内的零点数

$m_o$ ：单位圆外的零点数

$p_i$ ：单位圆内的极点数

$p_o$ ：单位圆外的极点数

$m_i+m_0=M;p_i+p_o=N$ ，在一个因果稳定的系统中，所有极点都位于单位圆内，即

$p_o=0,p_i=N$

现在我们令数字频率从0到 $2\pi$

$arg[\frac{H(e^{j\omega})}{K}]|_{\Delta\omega=2\pi}=2\pi(N-M)+2\pi(m_i-N)=2\pi(m_i-M)<0$

此时幅角变化为负，因此相位产生了延时，此时系统称为相位延时系统。相位角的变化大小与单位圆外的零点数目有关。

当零点全部位于单位圆内时，相位角变化最小， $m_i=M$

$arg[\frac{H(e^{j\omega})}{K}]|_{\Delta\omega=2\pi}=0$

此时成为最小相位延迟系统。

当零点全部位于单位圆外时，相位角变化最大， $m_o=M$

$arg[\frac{H(e^{j\omega})}{K}]|_{\Delta\omega=2\pi}=-2\pi M$

此时成为最大相位延迟系统。

对于延时系统的定义可以参看下表

FIR滤波器和IIR滤波器的延时

椭圆数字滤波器

求函数相关性

import numpy
import soundfile
from matplotlib import pyplot
from scipy import signal


def coherence(x, y, fs, n_fft, f_range):
    n_range = n_fft * f_range // fs
    n_move = n_fft // 2
    window = 0.5 * (1 - numpy.cos(2 * numpy.pi * numpy.arange(n_fft) / n_fft))
    n_frame = (len(x) - n_fft) // n_move + 1
    loc = 0
    p_xy, p_xx, p_yy = 0., 0., 0.
    for i in range(n_frame):
        fft_x = numpy.fft.rfft(window * x[loc:loc + n_fft])[:n_range]
        fft_y = numpy.fft.rfft(window * y[loc:loc + n_fft])[:n_range]
        p_xy = p_xy + fft_x.conj() * fft_y
        p_xx = p_xx + fft_x.conj() * fft_x
        p_yy = p_yy + fft_y.conj() * fft_y
        loc = loc + n_move
    c_xy = numpy.real(abs(p_xy) ** 2 / (p_yy * p_xx))
    f = numpy.fft.fftfreq(n_fft)[:n_range] * fs
    return f, c_xy



err, _ = soundfile.read("711-1.wav")
ref, sr = soundfile.read("ff-1.wav")

f2, coh = coherence(x1, y1, 16000, 4096, 20000)
att = 10 * numpy.log10(1 - coh)
fig, ax = pyplot.subplots(1, 1)
pyplot.xscale('log')
pyplot.xlim(10, 20000)
pyplot.plot(f2, att)
pyplot.show()

插值法

现在需要将两个系统响应相加，但是它们的横坐标不能对应，此时我们可以利用python的numpy工具包中的interp进行插值。

import numpy
import matplotlib.pyplot as plt

#先找到你要读入的txt文件路径
filepath_1 = 'filter.txt'
filepath_2=  '711 jianxiao.txt'


def loadmatrix(filename):
    returnMat = numpy.loadtxt(filename, delimiter='\t')  # 导入数据
    returnMat = returnMat[:, 0:3]  # 存放数据
    print(type(returnMat))
    return returnMat

def wrapto180(angle):
    return numpy.mod(angle + 180, 360) - 180


F_1=loadmatrix(filepath_1)
F_2=loadmatrix(filepath_2)

N=24#24分之一倍频程
f_i=20*numpy.exp2(numpy.arange(10*N)/N) #插值频率
f_p1=F_1[:,0]#f_primary 初始频率
SPL_p1=F_1[:,1]#SPL_primary 初始SPL
Phase_p1=F_1[:,2]#Phase_primary 初始相位
SPL_i1 = numpy.interp(f_i, f_p1, SPL_p1)  # 插值后SPL_interpolation
Phase_i1 = numpy.interp(f_i, f_p1, Phase_p1)  # 插值后Phase_interpolation

f_p2=F_2[:,0]#f_primary 初始频率
SPL_p2=F_2[:,1]#SPL_primary 初始SPL
Phase_p2=F_2[:,2]#Phase_primary 初始相位
SPL_i2 = numpy.interp(f_i, f_p2, SPL_p2)  # 插值后SPL_interpolation
Phase_i2 = numpy.interp(f_i, f_p2, Phase_p2)  # 插值后SPL_interpolation
SPL_i=SPL_i1+SPL_i2
Phase_i=Phase_i1+Phase_i2
# plt.plot(f_p1, SPL_p1, '-o')
plt.subplot(2,1,1)
plt.xlim(20,20000)
plt.semilogx(f_p1, SPL_p1, label='SPL_1')
plt.semilogx(f_p2, SPL_p2, label='SPL_2')
plt.semilogx(f_i, SPL_i,label='SPL_S')
# 添加图例
plt.legend()

plt.subplot(2,1,2)

plt.semilogx(f_p1, wrapto180(Phase_p1), label='Phase_1')
plt.semilogx(f_p2, wrapto180(Phase_p2), label='Phase_2')
plt.semilogx(f_i, wrapto180(Phase_i),label='Phase_S')
# plt.semilogx(f_p1, Phase_p1, label='Phase_1')
# plt.semilogx(f_p2, Phase_p2, label='Phase_2')
# plt.semilogx(f_i,Phase_i,label='Phase_S')
# 添加图例
plt.legend()
# 显示图像
plt.show(block=True)
plt.xlim(20,20000)

频谱的相关性计算

原文：

Python: scipy.signal.coherence的用法及代码示例_python coherence_嵌入式技术的博客-CSDN博客

信号处理中，信号的相关性是两种或两种以上信号分析的常用手段，用于寻找两种信号之间的关联性

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt

# 下面代码生成两个具有一共同特征的测试信号
fs = 16000				# 采样频率
N = 1e5					# 采样点数
amp = 50				# 信号幅值
freq = 1000				# 信号的特征频率
noise_power = 0.001 * fs / 2		# 信号的噪声功率
time = np.arange(N) / fs			# 信号的时间轴
b, a = signal.butter(2, 0.5, 'low')	# 构造 butter 滤波器
x = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
y = signal.lfilter(b, a, x)

x += amp*np.sin(2*np.pi*freq*time)
# y += np.random.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)
# plt.plot(time,x,color='b')
# plt.plot(time,y,color='r')
# plt.show()

# # 计算 x 与 y 的幅度平方相关估计
f, Cxy = signal.coherence(x, y, fs, nperseg=1024)

# 使用 welch 方法分别计算 x 与 y 各自的功率谱密度估计
f, Pxx = signal.welch(x, fs, nperseg=1024)
f, Pyy = signal.welch(y, fs, nperseg=1024)

# 下面绘制结果
lnxy = plt.semilogy(f, Cxy, color='r', label='$C_{xy}$')
plt.ylabel('Coherence')

# 绘制第二个 y 轴，用于标注 x 与 y 的功率谱密度估计
plt.twinx()
lnx = plt.semilogy(f, Pxx, color='g', label="$P_{xx}$")
lny = plt.semilogy(f, Pyy, color='b', label="$P_{yy}$")
plt.ylabel('PSD [$V^2/Hz$]')
plt.xlabel('frequency [Hz]')

# 绘制图例
lns = lnxy + lnx + lny
labs = [l.get_label() for l in lns]
plt.legend(lns, labs)
plt.show()

Butterworth滤波器使用

[b，a] = butter（n，Wn）返回归一化截止频率Wn的n阶低通数字巴特沃斯滤波器的传递函数系数。

b，a-传递函数系数

对于数字滤波器，传递函数可以表示为：

$\displaystyle H(z)=\frac{B(z)}{A(z)}=\frac{b(1)+b(2)z^{-1}+...+b(n+1)z^{-n}}{a(1)+a(2)z^{-1}+...+a(n+1)z^{-n}}$

对于模拟滤波器，传递函数可以表示为：

$\displaystyle H(s)=\frac{B(s)}{A(s)}=\frac{b(1)s^n+b(2)s^{n-1}+...+b(n+1)}{a(1)s^n+a(2)s^{n-1}+...+a(n+1)}$

[z，p，k] = butter（n，Wn，ftype）根据ftype的值和Wn的元素数来设计低通，高通，带通或带阻Butterworth滤波器，并返回其零，极点和增益。该语法可以包括先前语法中的任何输入参数。

对于数字滤波器，传递函数用z，p，k表示

$\displaystyle H(z)=k\frac{(1-z(1)z^{-1})(1-z(2)z^{-1})...(1-z(n)z^{-1})}{(1-p(1)z^{-1})(1-p(2)z^{-1})...(1-p(n)z^{-1})}$

对于模拟滤波器

$\displaystyle H(s)=k\frac{(s-z(1))(s-z(2))...(s-z(n))}{(s-p(1))(s-p(2))...(s-p(n))}$

使用python的butterworth滤波

from scipy.signal import butter, lfilter
import numpy
from matplotlib import pyplot as plt
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# Filter requirements.
order = 5
fs = 200      # sample rate, Hz
cutoff = 30  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)
t=numpy.arange(0,1,1/fs)
N=len(t)
# Demonstrate the use of the filter.
data = 1*numpy.sin(2*numpy.pi*20*t)+1*numpy.sin(2*numpy.pi*70*t)+numpy.random.normal(size=N)
y = butter_lowpass_filter(data, cutoff, fs, order)

# Frequency spectrum
H_1=numpy.fft.rfft(data)
H_2=numpy.fft.rfft(y)
f=numpy.fft.rfftfreq(len(t),1/fs)

# plot
plt.figure()

plt.subplot(2,1,1)
plt.plot(t,numpy.abs(data),'r-',label='Sin+Noise')
plt.plot(t,numpy.abs(y),'b--',label='Filtered')
plt.title('Time Spectrum')
plt.xlabel('Time [s]')
plt.ylabel('Magnitude')
plt.legend()

plt.subplot(2,1,2)
plt.plot(f,numpy.abs(H_1),'r-',label='Sin+Noise')
plt.plot(f,numpy.abs(H_2),'b--',label='Filtered')
plt.title('Frequency Spectrum')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude')
plt.legend()

plt.tight_layout()
plt.show()

电学端阻抗曲线

将力学端等效类比到电学端，画出阻抗曲线

import numpy as np
import matplotlib.pyplot as plt
Re=9.72 #电阻
Le=35e-4 #电感
# Le=35e-6
BL=0.483 #电力耦合系数
Rm=8e-3 #力阻
Mm=17.7e-6 #质量
Cm=23.65e-3 #力顺
f=np.arange(1,1000,5)
#
Z=Re+1j*2*np.pi*f*Le+(BL**2)/(Rm+1j*2*np.pi*f*Mm-1j*1/(2*np.pi*f*Cm))
plt.figure()
plt.plot(f,np.abs(Z))
plt.title('Impedance Curve')
plt.xlabel('Frequency [Hz]')
plt.ylabel('|Z|')
plt.show()