【笔记】高数学习笔记 - 区间再现 一、区间再现的表达$$ \int_a^bf(x)dx = \int_a^bf(a+b-x)dx = \frac{1}{2}\int_a^b[f(x)+f(a+b-x)]dx $$证明：$$ \int_a^bf(a+b-x)dx\overset{令t=a+b-x}{=}\int_b^af(t)(-1)dt=\int_a^bf(t)dt=\int_a^bf(x)dx $$例题$$ \begin{aligned} &1.1\qquad \int_0^\pi \frac{\sin x}{\sin x+\cos x}dx \\ &2.1\qquad \int_0^{\frac{\pi}{2}}\frac{1}{1+(\tan x)^\alpha}dx \\ &2.2\qquad \int_0^{+\inf}\frac{1}{(1+x^2)(1+x^{\alpha})}dx \\ &3.1\qquad \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}(\arctan e^x)\cdot \sin^2xdx \\ &4.1\qquad \int_{-2}^{2}x\cdot\ln(1+e^x)dx \\ &5.1\qquad \int_{\frac{\pi}{6}}^{\frac{\pi}{3}}\frac{\cos^2x}{x(\pi-2x)}dx \\ &6.0\qquad \int_0^{\frac{\pi}{4}}\frac{x}{\cos(\frac{\pi}{4}-x)\cdot\cos x}dx \\ &6.1\qquad \int\frac{1}{\sin(x+a)\cdot\sin(x+b)} \\ &6.2\qquad \int\frac{1}{\sin(x+3)\cdot\sin(x+5)} \\ &7.0\qquad \int_0^{\frac{\pi}{2}}\ln\sin xdx \\ &7.1\qquad\int_0^{\frac{\pi}{2}}\frac{x}{\tan x}dx \\ &7.2\qquad \int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{\cos x\ln \cos x}{1+\sin x+\cos x}dx \\ &8.0\qquad \int_0^1\frac{\ln(1+x)}{1+x^2}dx \\ &8.1\qquad \int_0^1\frac{\arctan x}{1+x}dx \\ &9.0\qquad \int_0^{n\pi}x\cdot |\sin x| dx \\ &9.1\qquad a_n=\int_0^{n\pi}x\cdot |\sin x| dx, 求\sum_{n=1}^{\inf}\frac{1}{a_n} \\ &10.0\qquad \int_0^1\frac{\arcsin\sqrt{x}}{\sqrt{x^2-x+1}}dx \\ &11.0\qquad \int_0^1x\cdot\arcsin 2\sqrt{x-x^2}dx \\ &12.0\qquad \int_0^1(1-x)^{100}\cdot xdx \\ &13.0\qquad \int_0^2x\cdot(x-1)\cdot(x-2)dx \\ &13.1\qquad \int_0^{2n}x\cdot(x-1)\cdot(x-2)\cdot\cdot\cdot\cdot(x-2n) dx \end{aligned} \\ $$例题解析之后再补吧。这是看凯哥的学习笔记。原视频可以参考：https://www.bilibili.com/video/BV1ah411Y789

【学习笔记】航天器在空中飞行的运动仿真 本文仅为某次课程大作业，有诸多不严谨、遗漏之处，全文仅供参考！一、大作业参考模型与数据某运载器从发射点Ls采用垂直发射方式发射，假定整个飞行过程中不存在级间分离，运载器与其有效载荷一同入轨或一同再入。假定运载器主动段只在其射面内飞行，不考虑侧向运动。假定地球为均质圆球且地球周围没有风。假定飞行过程中地心坐标系为惯性坐标系且其x轴指向平春分点。主动段纵向质心动力学方程组（书中式3.70）简化为如式（1）所示形式。$$ \begin{aligned} &v=\frac{P\cos(\alpha)-D}{m}-g\sin\theta \\ &\theta=\frac{P\cos(\alpha)+L}{mv}-g\frac{cos\theta}{v} \\ &x=v\cos\theta \\ &y=v\sin\theta \qquad\qquad\qquad\qquad\qquad\qquad\qquad(1)\\ &\phi=\theta+\alpha \\ &r=\sqrt{x^2+(y+R_e)^2} \\ &h=r-R_e \\ & m=m_0-mt \\ \end{aligned} $$其中 $\alpha$ 为飞行速度大小、$\theta$ 为速度倾角、$x,y$ 为运载器在发射坐标系下位置坐标、$\phi$ 为俯仰角、$r$ 为地心距、$h$ 为飞行高度、$P$ 为总的发动机推力大小、$\alpha$ 为攻角、$g$ 为引力加速度大小、$L$ 为气动升力大小、$D$ 为气动阻力大小、$R_e$ 为地球半径、$m$ 为运载器当前质量、$m_0$ 为运载器起飞质量、$m$ 为质量秒耗量、$t$ 为飞行时间。a)大气模型：大气密度采用指数模型，大气压强采用拟合模型如式 $(2)$ 所示，假定声速 $v_s$ 不随飞行高度变化，$v_s=340m/s$ 为常值。$$ p(h)= \begin{cases} p_0e^{-0.142h}, & h\le h\le 92km \quad 其中p_0=101325Pa\\ p(h_{pf}), & h>92km \quad 其中h_{pf}=92km \end{cases} \tag{2} $$b)引力模型：采用平方反比模型，如式 $(3)$ 所示。为地球引力常数。$$ g=\frac{\mu}{r^2}\tag{3} $$c)气动模型：采用拟合模型如下所示。（$M$ 为马赫数、$C_L$ 为升力系数、$C_D$ 为阻力系数、注意此拟合模型中攻角 $\alpha$ 单位为度）$$ \begin{aligned} &C_L(\alpha,M)=-0.0005225\alpha^2+0.03506\alpha-0.04857M+0.1577 \\ &C_D(\alpha,M)=0.0001432\alpha^2+0.00558\alpha-0.01048M+0.2204 \end{aligned} $$d)主动段飞行程序角分段设计，如 $(4)$ 式所示。$$ \phi(h)= \begin{cases} \frac{\pi}{2} &t\le t_1 \\ a(t-t_2)^2+\phi_f \quad 其中a=\frac{\frac{\pi}{2}-\phi_f}{(t_1-t_2)^2} & t_1 < t \le t_2 \\ \phi_f & t>t_2 \end{cases} \tag{4} $$其中，$t=\sqrt{40/(\frac{P_{0总}}{m_0*g_0}-1)}$, $P_{0总}$ 为总的地面推力，地面引力加速度 $g_0=9.8 m/s^2$。$\phi_f$、$t_2$ 为给定设计参数（如下）。e)仿真计算采用数据地球半径 $R_E=6378.1km$、飞行器参考面积（特征面积）为 $20m^2$、起飞质量 $m_0=200000kg$ 、飞行器空重即结构质量为 $90000kg$。主发动机共 $3$ 台，单台发动机地面推力为为 $1500000N$、发动机喷口截面积为 $0.05 m^2$。主动段飞行采用耗尽关机方式（即将所有燃料全部用完发动机才关机）。发射点Ls经纬度分别为 $\lambda_0=60\deg, \phi_0=45\deg$，地心方位角 $\alpha_0=30\deg$。假定在发射瞬时 $t=0$ 时刻，飞行速度 $v_0=0.1m/s$、速度倾角 $\theta_0=90\deg$、运载器在发射坐标系下位置坐标 $x_0=0m, y_0=0m$。二、大作业任务第一问题目按照上述模型与数据，分别在：$$ \begin{aligned} &Case1: \quad m=500kg/s，t2=200s，\phi_f=20\deg（换成弧度） \\ &Case2: \quad m=300kg/s，t2=230s，\phi_f=5\deg（换成弧度） \end{aligned} $$两种情况下通过编写主动段飞行仿真程序（数值积分式（1）），进行主动段轨迹仿真。给出仿真结果（飞行状态如位置、速度、高度、欧拉角等随时间变化曲线），得出两种情况下主动段终点 $k$ 处在发射坐标系下的飞行速度、速度倾角、位置坐标。解答由题意，弹道的主动段飞行轨迹式一阶变系数方程组，可以用数值的方法求得数值解。本题可以通过欧拉法的方式进行求解。求解可以使用随后列出的仿真程序。在程序中输入所需的，t2，的值以及合适的步长，即可求出所需的数值解。以下为在设定不同的步长时仿真所的出的case1的x坐标值、所花费的时间与步长的关系图：观察可发现，当步长设定为0.001时计算花费了。。的时间，且在此之后所得出的数值差距不大。因此之后的仿真均采用0.001作为步长进行计算。经过计算，可得在case1与case2下主动段终点k处的飞行速度、飞行倾角、位置坐标分别为：并绘出飞行状态随时间的变化曲线（对于每个case的曲线分别为：速度-时间变化曲线、飞行高度-时间变化曲线、y-x曲线、攻角-时间变化曲线、速度倾角-时间变化曲线、俯仰角-时间变化曲线）：附：计算所需的程序：# homework.py import numpy as np import matplotlib.pyplot as plt # basic constants PI = 3.14159265 beginMass = 200000 # mkg beginGravity = 9.8 structureMass = 90000 Sref = 20 Propulsion = 1500000 * 3 # m3 * 1.5e6 N baseDensity = 1.125 # mkg/m^3 baseAcousitc = 340 # mm/s Rearth = 6378100 # m mu = 398600000000000 # m^3/s^2 def getGravity(radius): return mu / radius / radius def getAtmosphericDensity(height): return baseDensity * np.exp(-1.406 * 10 ** (-4) * height) def getAerodynamicDrag(height, alpha, velocity): density = getAtmosphericDensity(height) mach = velocity / baseAcousitc alpha_deg = alpha * 180 / 3.14159265 coefficient = 0.0001432 * alpha_deg ** 2 + 0.00558 * alpha_deg - 0.01048 * mach + 0.2204 return 0.5 * coefficient * density * velocity **2 * Sref def getAerodynamicLift(height, alpha, velocity): density = getAtmosphericDensity(height) mach = velocity / baseAcousitc alpha_deg = alpha * 180 / 3.14159265 coefficient = -0.0005225 * alpha_deg ** 2 + 0.03506 * alpha_deg - 0.04857 * mach + 0.1577 return 0.5 * coefficient * density * velocity **2 * Sref def getPhi(time, t1, t2, phiF): if(time < t1): return PI / 2 elif (time < t2): a = (PI / 2 - phiF) / (t1-t2)**2 return a * (time-t2)**2 + phiF else : return phiF def run(dm, t2, phiF, stepLength): # calculate constant t1 = np.sqrt(40/(Propulsion / (beginMass * beginGravity) - 1)) # calculate variables nowTime = 0.0 mass = beginMass v = 0.1 x = 0 y = 0 theta = PI / 2 dataSet = [] cnt = 0 while(True): lengthToEarthCenter = np.sqrt(x ** 2 + (y + Rearth) ** 2) height = lengthToEarthCenter - Rearth phi = getPhi(nowTime, t1, t2, phiF) alpha = phi - theta dv = (Propulsion * np.cos(alpha) - getAerodynamicDrag(height, alpha, v)) / mass - getGravity(lengthToEarthCenter) * np.sin(theta) dtheta = (Propulsion * np.sin(alpha) + getAerodynamicLift(height, alpha, v)) / (mass * v) - getGravity(lengthToEarthCenter) * np.cos(theta) / v v += dv * stepLength theta += dtheta * stepLength x += v * np.cos(theta) * stepLength y += v * np.sin(theta) * stepLength mass -= dm * stepLength nowTime += stepLength cnt += 1 if(mass <= structureMass): break if(cnt % 1000 == 0): dataSet.append([nowTime, alpha, theta, phi, x, y, height, v]) print('计算结束！') print('位置(xk, yk, zk) = ('+str(x/1000)+'km, '+str(y/1000)+'km, 0km)、速度vk='+str(v)+'、速度倾角θk='+str((theta+np.arctan(x/(y+Rearth)))*180/3.14159265)+'°') return dataSet def showData(datas): for data in datas: time = [el[0] for el in data] alpha = [el[1] for el in data] theta = [el[2] for el in data] phi = [el[3] for el in data] x = [el[4] for el in data] y = [el[5] for el in data] height = [el[6] for el in data] v = [el[7] for el in data] # 使用子图 plt.figure() plt.suptitle('The plot of case '+str(cnt), size=16) # 整张图的总标题 plt.subplot(2, 3, 1) # （行，列，活跃区） plt.title('velocity - time') # 子图的小标题 plt.plot(time, v, ls="-", lw=2) plt.subplot(2, 3, 2) plt.title('height - time') # 子图的小标题 plt.plot(time, height, ls="-", lw=2) plt.subplot(2, 3, 3) plt.title('y - x') # 子图的小标题 plt.plot(x, y, ls="-", lw=2) plt.subplot(2, 3, 4) plt.title('alpha - time') # 子图的小标题 plt.plot(time, alpha, ls="-", lw=2) plt.subplot(2, 3, 5) plt.title('theta - time') # 子图的小标题 plt.plot(time, theta, ls="-", lw=2) plt.subplot(2, 3, 6) plt.title('phi - time') # 子图的小标题 plt.plot(time, phi, ls="-", lw=2) plt.show() print('case1:') case1 = run(500, 200, 20.0 / 180.0 * 3.14159265, 10**(0)) print('case2:') case2 = run(300, 230, 5.0 / 180.0 * 3.14159265, 0.35) showData([case1, case2])程序的输出为：case1:位置(xk, yk, zk) = (403.4326437092588km, 252.27295090009528km, 0km)、速度vk=5698.847570444022、速度倾角θk=23.51929491334418°case2:位置(xk, yk, zk) = (1449.6315824637384km, 228.72349043694848km, 0km)、速度vk=10529.850447985076、速度倾角θk=12.21735847248481°第二问case1：根据（1）计算结果，有：$$ \begin{aligned} & e=\sqrt{1+\nu_k(\nu_k-2)\cos^2\Theta_k}=0.8664258403313946 \\ & p=r^2_kv^2_k\cos^2\Theta_k = r_k\nu_kcos^2\Theta_k=12014996 \end{aligned} $$再因为$$ \begin{aligned} & r_k=6722273.510508056\, m\\ & r_p=6437435.844933111\, m \\ & r_l=6428100\, m \\ \end{aligned} $$因此有：$$ r_k\geq r_p\geq r_l \quad成立 \qquad(1) $$再将$\Theta_k$，$r_l$，$r_k$，$v_k$带入，可得：$$ \begin{aligned} & \cos\Theta_k=0.916925739930821\\ & \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}=1.0256422532299725 \\ & v^2_k=32476863.63115573 \\ & \frac{2\mu r_l}{1+r_k/r_l}=59021397.478749804 \\ \\ & \cos\Theta_k \geq \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}\quad不成立\qquad(2) \\ & v^2_k\geq \frac{2\mu r_l}{1+r_k/r_l}\quad不成立\qquad(3) \end{aligned} $$因为(2), (3)不成立，因此此情况下能成为导弹。case2：根据（1）计算结果，有：$$ \begin{aligned} & e=\sqrt{1+\nu_k(\nu_k-2)\cos^2\Theta_k}=0.8664258403313946 \\ & p=r^2_kv^2_k\cos^2\Theta_k = r_k\nu_kcos^2\Theta_k=12014996.606458724 \end{aligned} $$将$x_k$，$y_t$，$R_e$代入，可得：$$ r_k=6722273.510508056 \\ r_p=6437435.844933111 \\ r_l=6428100 \\ \quad\\ r_k\geq r_p\geq r_l \quad成立 \qquad(1) $$再将$\Theta_k$，$r_l$，$r_k$，$v_k$带入，可得：$$ \begin{aligned} & \cos\Theta_k=0.9773518259937374\\ & \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}=0.9763775707871245 \\ & v^2_k=110877750.45693152 \\ & \frac{2\mu r_l}{1+r_k/r_l}=57429289.15560445 \\ \\ & \cos\Theta_k \geq \frac{r_l}{r_k}\sqrt{1+\frac{2\mu}{v^2_k}(\frac{1}{r_l}-\frac{1}{r_k})}\quad成立\qquad(2) \\ & v^2_k\geq \frac{2\mu r_l}{1+r_k/r_l}\quad成立\qquad(3) \end{aligned} $$因为(1), (2), (3)均成立，因此此情况下能成为卫星。(3).解：由（2）可知，case1的情况可以成为导弹。此时，根据（1）可得主动段终点的参数：$$ \begin{aligned} & \Theta_k=23.51929491334418 \\ & \nu_k=0.5412241844023746 \\ & r_k=6642635.272693954 \\ \end{aligned} $$代入被动段射程公式，可得：$$ \begin{aligned} & A=2R(1+\tan^2\Theta_k)-\nu_k(R+r_k)=8125218.006125015 \\ & B=2\nu_k R\tan\Theta_k=3004693.9006598108 \\ & C=\nu _k(R-r_k)=-143172.88720944524 \end{aligned} $$因此有：$$ A\tan^2\frac{\beta_c}{2}-B\tan\frac{\beta_c}{2}+C=0 $$再因为$\beta\geq0$，因此有：$$ \begin{aligned} & \tan\frac{\beta_C}{2}=\frac{B+\sqrt{B^2-4AD}}{2A}=0.36979822647587485 \\ & \quad \\ & L_{KC}=R\beta_C=R\arctan(0.7084048610762707)=4518277.044430562 \end{aligned} $$同理，对自由段，有：$$ \begin{aligned} & A=2r_k(1+\tan^2\Theta_k)-2r_k\nu_k=8611327.040417206 \\ & B=2r_k\nu_k \tan\Theta_k=3129315.264447287 \\ & C=0 \\ & \quad \\ & \tan\frac{\beta_E}{2}=\frac{B}{A}=\frac{\nu_k\sin\Theta_k\cos\Theta_k}{1-\nu_k\cos^2\Theta_k} \\ & \quad \\ & L_{KE}=R\beta_E=R\arctan(0.6971160493800195)=4446275.874550702 \end{aligned} $$综上所述被动段的射程$L_{KC}=4518277.044430562$，自由段的射程$L_{KE}=4446275.874550702$。$$ \begin{aligned} & \tan\Theta_{K·OPT}=\sqrt{\frac{\nu_k[2R-\nu_k(R+r_k)]}{2R\nu_k-4(R-r_k)}}=0.6229555659247278 \\ & \tan\frac{\beta_{Cmax}}{2}=\sqrt{\frac{\nu_k[R\nu_k-2(R-r_k)]}{2[2R-\nu_k(R+r_k)]}}=0.4344003120021654 \\ & \quad \\ & \tan\Theta_{KE·OPT}=\sqrt{1-\nu_k}=0.6773299163610192 \\ & \tan\frac{\beta_{Emax}}{2}=\frac{1}{2}\frac{\nu_k}{\sqrt{1-\nu_k}}=0.3995277421895981 \end{aligned} $$所以，主动段终点k处的当地速度倾角不是最优速度倾角。（4）解：由题意可得导弹的最小负加速度：$$ \begin{aligned} & h_m=\frac{1}{\beta}\ln(-\frac{C_xS_m\rho_0}{m\beta \sin\Theta_e})=9040.165079085044 \\ & v_m=v_ee^{-\frac{1}{2}}=3456.5257765031506 \\ & \dot{v_m}\frac{\beta v^2_e}{2e}\sin\Theta_e=1571.31518785745 \end{aligned} $$导弹的最大驻点热流：$$ \begin{aligned} & h_{m2}=\frac{1}{\beta}\ln(-\frac{3C_xS_m\rho_0}{m\beta \sin\Theta_e})=16853.908241731628 \\ & v_{m2}=v_ee^{-\frac{1}{6}}=4823.970321318141 \\ & (q_s)_{max}=k_s\sqrt{\frac{-m\beta\sin\Theta_e}{3eCxSm}}v^3_e=11825312565572.254 \end{aligned} $$导弹的再入过程总吸热量：$$ \begin{aligned} & v_c=v_ee^{\frac{2B}{\beta\sin\Theta_e}\rho_e}=161.3273896537413 \\ & Q=\frac{mc_f^{\prime}S_T}{4C_xS_M}(v^2_e-v^2_c)=9126797935641.537 \end{aligned} $$分析最大驻点热流$(q_s)_{max}$可以发现，若要降低最大驻点热流，应减小再入角$|\Theta_e|$，但如果$|\Theta_e|$过小，又会增加飞行的时间，使总吸热量增加。因此，调整时需要兼顾总吸热量。合理的选择一个再入角$\Theta_e$，是弹道再入的一个重要问题。（5）解：由（2）可知，case2的情况可以成为卫星。此时，根据（1）可得主动段终点的参数：由上述参数可求出弹道关机点时的卫星地心矩矢量r速度矢量v分别为：$$ \begin{aligned} & r=[-369797.99601624 \quad 296122.06928999 \quad -44279.81467744]^T \\ & v=[-4513.92140207 \quad 3078.23939085 \quad -1620.43804561]^T \end{aligned} $$并求得其矢积：$$ \begin{aligned} h&=r\times v \\ &=[-3.43543597e+08 \quad -3.99359139e+08 \quad 1.98344988e+08] \end{aligned} $$由卫星轨道与赤道惯性坐标系之间的关系可得动量矩h的三个分量与$i$，$\Omega$的关系：$$ \begin{aligned} & h_x=h\sin i\sin \Omega \\ & h_y=-h\sin i\cos \Omega \\ & h_z=h\cos i \end{aligned} $$因此易得：$$ \begin{aligned} & i=\arccos(\frac{h_z}{h})=1.2106983360402104\,rad \\ & \Omega=\arctan(\frac{-h_z}{h_y})=0.46097059130336926\,rad \end{aligned} $$由开普勒定律可知，卫星绕地球的运动轨迹为一椭圆。因此易得：$$ \begin{aligned} & e=0.8664258403313946 \\ & a=\frac{h^2}{\mu(1-e^2)}= \end{aligned} $$再因为卫星轨道方程的标准极坐标形式：$$ \begin{aligned} & r=\frac{p}{1+e\cos(\theta-\omega)} \\ & p=h^2/\mu \end{aligned} $$将关机点参数$r=, \theta=0$代入，可解得近地点幅角：$$ \omega= $$由几何关系可得出卫星轨道的偏近点角$E$与真近点角$f$之间的关系：$$ a\cos E=ae+r\cos f \\ b\sin E=r\sin f $$已知卫星的位置与时间关系的开普勒方程：$$ t-t_p=\sqrt{\frac{a^3}{\mu}}(E-e\sin E) $$将关机点参数$f=-\omega, t=0$代入，可解得卫星经过近地点的时刻$t_p$：$$ t_p=-\sqrt{\frac{a^3}{\mu}}(E-e\sin E)= $$至此，给出了卫星运动方程的6个积分常数。References忘了（

【学习笔记】FourierTransform-关于二维DCT变换可以被表示为矩阵相乘这档事 前言暂略。前置知识需要一定编程基础，理解连续/离散傅里叶变换原理与计算方法。正文首先，让我们来看一下DCT的公式：$$ F(u)= \begin{cases} \frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}f(x), & n\ne m \\ \frac{2}{\sqrt{N}}\sum_{x=0}^{N-1}f(x)\cos\frac{(2x+1)u\pi}{2N}, &n=m \end{cases} $$具体推导过程可见前文。上式即为DCT的公式表达。而在具体到计算时，例如要处理一个长度为8的时域信号时，可以展开为：通过观察我们可以发现在变换的具体表达上很像矩阵相乘的样子。因此，我们可以进一步改写为：其中，最左边的$[F(u)]$为变换后的矩阵，即为频域的矩阵。最右边的$f[x]$为变换前的时域矩阵。而中间的矩阵在N确定时为一确定矩阵，我们把它称之为DCT的变换矩阵A。而A可以写为： $$ A[x,y]= \begin{cases} \frac{1}{\sqrt{N}}, & y= 0 \\ \frac{2}{\sqrt{N}}\cos\frac{(2x+1)y\pi}{2N}, &y\ne 0 \end{cases} $$所以，我们可以把一维的DCT变换表现为矩阵相乘的形式，即：$F=A[f(x)]$。这是一维的情况，那么二维的呢？我们首先来看下二维离散余弦变换的一般公式（为了写的方便稍微换了下写法）：$$ F(u,v)=\alpha(u)\alpha(v)\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)\cos\frac{(2x+1)u\pi}{2M}\cos\frac{(2y+1)v\pi}{2N} $$其中：$$ \alpha(u)= \begin{cases} \frac{1}{\sqrt{N}}, & u= 0 \\ \frac{2}{\sqrt{N}}, &u\ne 0 \end{cases} \alpha(v)= \begin{cases} \frac{1}{\sqrt{N}}, & v= 0 \\ \frac{2}{\sqrt{N}}, &v\ne 0 \end{cases} $$可以看到，二维的DCT表达起来和一维的DCT变换差不多，基本上就是在原始数据的x方向做一次一维DCT，然后再在y方向再做一次一维 DCT。而事实上也就是那样。所以对于一个指定的$[f(x,y)]$，我们要做的就是对它的每一列都做一次DCT变换，再转置一次后再做一次DCT变换，最后再转置回来，就完成了二维的DCT。具体可以表示为：$$ [F(u,v)]=(A((A[f(x,y)])^T)) ^T $$再因为转置矩阵的性质：$$ (AB)^T = B^T A^T $$所以我们可以进一步化简：$$ [F(u,v)]=(A((A[f(x,y)])^T)) ^T\\ =(A([f(x,y)]^T A^T)) ^T\\ =([f(x,y)]^T A^T) ^T A^T\\ =A[f(x,y)]A^T $$(写的格式有点问题见谅。。）所以我们就得到了二维DCT变换公式：$$ F(u,v)=A[f(x,y)]A^T $$编码实现为了进一步巩固学习内容，这里使用C语言写一个小DEMO，实现二维DCT的效果。代码除了math库和stdio库外不包含任何第三方库。编写完成后会使用matlab对所得结果进行验证。首先，进行DCT需要用到矩阵运算，这里先给出简易的矩阵运算代码。typedef struct { int rows; int cols; double** p; } Matrix; Matrix initMatrix(int rows, int cols) { Matrix mat; mat.rows = rows; mat.cols = cols; mat.p = (double**)malloc(sizeof(double*)*rows); // new double*[rows];//分配rows_num个指针 for (int i = 0; i < rows; ++i) { mat.p[i] = (double*)malloc(sizeof(double)*cols); // new double[cols];//为p[i]进行动态内存分配，大小为cols } return mat; } Matrix initMatrixWithValue(int rows, int cols, double* p) { Matrix mat = initMatrix(rows, cols); for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { mat.p[i][j] = *p++; } } return mat; } Matrix matrixTrans(Matrix mat) { Matrix res = initMatrix(mat.cols, mat.rows); for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { res.p[j][i] = mat.p[i][j]; } } return res; } Matrix matrixMul(Matrix m1, Matrix m2){ Matrix res = initMatrix(m1.rows, m2.cols); for (int i = 0; i < m1.rows; i++) { for (int j = 0; j < m2.cols; j++) { res.p[i][j] = 0; for (int k = 0; k < m1.cols; k++) { res.p[i][j] += (m1.p[i][k] * m2.p[k][j]); } } } return res; } //矩阵显示 void showMatrix(Matrix mat) { //cout << rows_num <<" "<<cols_num<< endl;//显示矩阵的行数和列数 for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { printf("%.4f ", mat.p[i][j]); } printf("\n"); } printf("\n"); }上述代码实现了矩阵的初始化，加法，乘法和转置运算。Matrix get4x4DCTMatrixL(){ Matrix mat = initMatrix(4, 4); for(int i = 0;i < 4;i++){ for(int j = 0;j < 4;j++){ double c = (i==0) ? sqrt(1.0/4) : sqrt(2.0/4); mat.p[i][j] = c * cos( (j + 0.5) * PI * i / 4); } } return mat; }这段代码即为生成DCT变换矩阵的代码。这里是一个处理4x4的数据的矩阵。之后若是要进行DCT变换，只要按照上文去做计算即可。为了进行测试，这里我先用matlab生成一段数据用于测试。代码如下：clc;clear f = (rand(4,4)*100); dct2(f);f = 40.1808 18.3908 90.2716 33.7719 7.5967 23.9953 94.4787 90.0054 23.9916 41.7267 49.0864 36.9247 12.3319 4.9654 48.9253 11.1203 ans = 156.9409 -54.8586 -28.9792 51.3964 43.0922 -19.6215 -0.1741 18.9064 -26.9619 28.4906 -3.5949 26.3422 -6.7750 39.6837 -3.5258 -9.3417接下来我们就在c语言中对上述数据进行计算，并对答案进行比较。代码如下：#include <stdio.h> #include <stdlib.h> #include <math.h> #define PI 3.1415926f typedef struct { int rows; int cols; double** p; } Matrix; Matrix initMatrix(int rows, int cols) { Matrix mat; mat.rows = rows; mat.cols = cols; mat.p = (double**)malloc(sizeof(double*)*rows); // new double*[rows];//分配rows_num个指针 for (int i = 0; i < rows; ++i) { mat.p[i] = (double*)malloc(sizeof(double)*cols); // new double[cols];//为p[i]进行动态内存分配，大小为cols } return mat; } Matrix initMatrixWithValue(int rows, int cols, double* p) { Matrix mat = initMatrix(rows, cols); for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { mat.p[i][j] = *p++; } } return mat; } Matrix matrixTrans(Matrix mat) { Matrix res = initMatrix(mat.cols, mat.rows); for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { res.p[j][i] = mat.p[i][j]; } } return res; } Matrix matrixMul(Matrix m1, Matrix m2){ Matrix res = initMatrix(m1.rows, m2.cols); for (int i = 0; i < m1.rows; i++) { for (int j = 0; j < m2.cols; j++) { res.p[i][j] = 0; for (int k = 0; k < m1.cols; k++) { res.p[i][j] += (m1.p[i][k] * m2.p[k][j]); } } } return res; } //矩阵显示 void showMatrix(Matrix mat) { //cout << rows_num <<" "<<cols_num<< endl;//显示矩阵的行数和列数 for (int i = 0; i < mat.rows; i++) { for (int j = 0; j < mat.cols; j++) { printf("%.4f ", mat.p[i][j]); } printf("\n"); } printf("\n"); } Matrix get4x4DCTMatrixL(){ Matrix mat = initMatrix(4, 4); for(int i = 0;i < 4;i++){ for(int j = 0;j < 4;j++){ double c = (i==0) ? sqrt(1.0/4) : sqrt(2.0/4); mat.p[i][j] = c * cos( (j + 0.5) * PI * i / 4); } } return mat; } int main(){ Matrix dctm = get4x4DCTMatrixL(); double testF_a[4*4] = { 40.1808, 18.3908, 90.2716, 33.7719, 7.5967, 23.9953, 94.4787, 90.0054, 23.9916, 41.7267, 49.0864, 36.9247, 12.3319, 4.9654, 48.9253, 11.1203 }; Matrix testF = initMatrixWithValue(4, 4, testF_a); Matrix tmp = matrixMul(dctm, testF); Matrix ans = matrixMul(tmp, matrixTrans(dctm)); showMatrix(ans); return 0; }我们把它编译运行试试：[ame@RaMizy ws]$ gcc test_dct.c -o test_dct -lm && ./test_dct 156.9409 -54.8586 -28.9792 51.3964 43.0922 -19.6215 -0.1741 18.9063 -26.9619 28.4906 -3.5949 26.3423 -6.7750 39.6837 -3.5258 -9.3417 可以看到，我们的代码计算结果与matlab的结果在合理的误差范围内基本一致。结果存在的微小误差可认为来自于不同的计算方式及精度带来的截断误差。结尾至此，已经基本讲完了DCT的矩阵相关。但由于时间所限，我也不可能讲的很完整，也难免会有些问题。就之后再补充吧。原本也想讲下像是蝶式算法之类的快速算法，就之后再细写吧。reference先放这里了。Reference二维离散余弦变换（2D-DCT）DCT变换展开时用到的示例图片DCT变换自学笔记离散余弦变换（DCT）的来龙去脉图像的DCT算法 matlab代码参考快速离散余弦变换代码实现(FDCT)H.264中整数DCT变换，量化，反量化，反DCT究竟是如何实现的？H264___DCT蝶形算法____理解H264-整数DCT变换和蝶形变换代码实现