MATLAB编写m函数理解 y=f(g(x))*h(x)

qq_45111955

已于 2022-10-23 09:19:07 修改

阅读量589

点赞数

分类专栏： matlab语法文章标签： matlab 算法

于 2022-10-22 17:07:20 首次发布

本文链接：https://blog.csdn.net/qq_45111955/article/details/127463260

版权

matlab语法专栏收录该内容

21 篇文章 3 订阅

订阅专栏

2.1 简单函数作用矩阵=矩阵，代替for循环

任务：

将数学形式： $y=f(g(x))$ ，编写为MATLAB m函数。

注：

1.任何字母本质上自变量都是k

MATLAB中，任何字母本质上自变量都是k，k表示第k个数。即任何函数再复杂，再复合，都可以for循环笨方法一个一个算出来。

$y=[y(k)]=[f(g( x(k) ))]$ ，y,x,[]表示矩阵，y(k),x(k)表示一个数。

则MATLAB里，

k=1:K

x=[x(k)]=x1:xs:x2

$y(k)=f(g(x(k)))$

+for 循环

end

但下标自变量k只能从1开始，不能从0。但字母值可以从零开始，如n=0:10;

2.除了笨方法，可能可以用特别的技巧：

2.1 简单函数作用矩阵=矩阵，代替for循环

若数学形式：y=f(x)，f都是简单的函数，比如exp,sin等，则MATLAB里可：

y=f(x) % y,x均为矩阵

不用

for k..

y(k)=f(x(k))

end

2.2 矩阵点乘替代函数相乘

若数学形式：y=f(x)g(x)，f,g都是简单的函数，比如exp,sin等，则MATLAB里或许可以省去for循环，通过：

y=f(x) .* g(x) % y,x均为矩阵。

2.3 矩阵相乘替代求和符号

若数学形式： $y=\sum_{k=1}^{N}f(x(k)))\cdot g((x(k)))$ ，则MATLAB里或许可以矩阵相乘替代求和符号，通过：

f=[f(k)]=[f(x(k))]

g=[g(k)]=[g(x(k))]

y=f * g'，f和g是k所控制的最高层变量。

3.例题：DTFT *Ts

求离散时间傅里叶变换 X(exp(jw))* Ts

编程出： $X(e^{j\Omega T_{s}})=\sum_{n=0}^{N-1}x(n)e^{-j\Omega T_{s}n} \cdot T_{s}$ 。

解：

1. X自变量看似exp(j*omega*T)，很复杂，但其实是X=[X(k)]。

采用逻辑自底层而上。

写成好看的形式：

$\Omega ({\color{Red} k})=\Omega1:\Delta :\Omega2$

$x=[x({\color{Red} k})]=e^{-j\Omega T_{s} n({\color{Red} k})}cos(w0 * T*n({\color{Red} k})))$ 。

$X({\color{Red} k})=\sum_{n=0}^{N-1}x(n({\color{Blue} k'})))e^{-j\Omega({\color{Red} k}) T_{s}n({\color{Blue} k'})} \cdot T_{s}$ 其实k'没必要标了，只要写到自变量（k），即

$X({\color{Red} k})=\sum_{n=0}^{N-1}x([n])e^{-j\Omega({\color{Red} k}) T_{s}[n]} \cdot T_{s}$ ,

2. exp,cos为简单函数，MATLAB里可以省去for，直接作用矩阵产生矩阵。

3. exp * cos 简单函数相乘，MATLAB里可以点乘 .* 。

$\begin{array}{c} x = [x(k)] =[{e^{ - j\Omega {T_s}n({\color{Red} k})}}cos(w0*T*n({\color{Red} k})))]\\ {\rm{ = }}{e^{ - j\Omega {T_s}{\color{Green} [n]}}}.*cos(w0*T*{\color{Green} [n]})) \end{array}\:$

4. 有求和符号，MATLAB可以通过矩阵相乘简化。

$\begin{array}{c} X({\color{Red} k}) = \sum\limits_{n = 0}^{N - 1} x (n({\color{Blue} k'}))){e^{ - j\Omega ({\color{Red} k}){T_s}n({\color{Blue} k'})}} \cdot {T_s}\\ {\rm{ = }}{T_s}{e^{ - j\Omega ({\color{Red} k}){T_s}[n]}}*[x]' \end{array}$

代码：

% 复合函数
% 求DTFT*Ts
clc;clear;close all;


%% x
% 参数
fs=3;   %采样频率
f0=0.5;
omega0=2*pi*f0;
Ts=1/fs;
T=100; %采样时间
N=T*fs;
%自变量
n=0:N-1;
zeta=0.05;
x=exp(-zeta*Ts*[n]) .* cos(omega0*Ts*[n]);

%% X
omega=[0: 0.01: 4*pi] /Ts;    %omega不必与[x(k)]数目相同，因为DTFD的w是连续的
for k=1:length(omega)

    X(k)=Ts * exp(-1i*omega(k)*Ts*[n]) * [x]';

end

% 连续时间傅里叶变换 FT
Y=1/2*(1./(zeta+1i*(omega-  omega0))+1./(zeta+1i*(omega+omega0)));
 
h=figure(1);
set(h,'Position',[150,150,800,500])
semilogy(omega/2/pi,abs(Y),'b','LineWidth',2)
xlabel('\Omega/2 \pi =f','FontName','Times New /*不要回车换行*/ Roman','FontSize',14)
ylabel('Amplitude','FontName','Times New /*不要回车换行*/ Roman','FontSize',14)
grid on

hold on 
semilogy(omega/2/pi,abs(X),'k','LineWidth',2)
set(gca,'FontName','Times New Roman','FontSize',14)
legend('True,FT',['DTFT*Ts,f_s=',num2str(fs)])

仿真结果：