Matlab用法记录

1、把数据写到txt文件

fp = fopen('filename.txt', 'w'); // 打开名为filename的txt文件，如果文件不存在就自动创建
fprintf(fp, '%d\n', data); //把整型变量data写入txt文件
fclose(fp);

2、保留小数

floor(x)：向下取整
ceil(x): 向上取整
round(x)：取最接近的整数
fix(x)：直接取整数部分

3、plot3

  下面代码相当于绘制点(1,2,3)、(2,3,4)和(3,4,5)。所以x，y和z必须是相同长度的向量。

x=1:3;
y=2:4;
z=3:5;
plot3(x,y,z);

4、meshgrid

  meshgrid()能够生成网格。网格怎么理解？以下面代码为例。X和Y为矩阵，X的每一行都是x，Y的每一列都是y，那么就得到了一个2x3网格。这个网格上有点(1,3)、(1,4)、(1,5)、(2,3)、(2,4)和(2,5)。那么Z必须是一个2x3的矩阵，(i,j,Zij)为三维空间中的点。

x=1:2;
y=3:5;
[X,Y]=meshgrid(x,y);
Z=X+Y;
surf(X,Y,Z);

5、surf

6、打印到命令行

  第一种是disp(变量名)，第二种是fprintf(格式化字符串, …）

7、else if

  Matlab用的是elseif，合起来的。

8、实现do…while

while(expression1)
	statement
	if(expression2)
		break;
	end
end

9、求两个向量的距离

dis=sqrt(nomr(a-b));

10、产生随机数（噪声）

  法一：unifrnd()

unifrnd(a,b,m,n); % mxn阶的[a,b]均匀分布
unifrnd(a,b); % [a,b]均匀分布

  法二：rand()（产生均匀噪声）

rand(m,n); % mxn阶的[0,1]均匀分布
rand(n); % nxn阶的[0,1]均匀分布

  法三：randn()（产生高斯噪声）

randn(m,n);	% mxn阶的高斯分布（数的大小为-∞到＋∞，大概率落在-3到3之间，0的概率最大）
randn(n);	% nxn阶的高斯分布

11、符号函数

syms x
f(x)=x^2;
fplot(f(x));
y=f(2);

12、函数

function [a,b] = fun_name(x1,x2)

  记录输入变量个数，用nargin

if(nargin==1)
	% do something
elseif(nargin==2)
	% do something
else
	% do something
end

13、全局变量

  在所有需要用到全局变量的.m文件里都要声明一下。

function func()
global x1
...

14、写函数

%采样
T=1;                %s        %观测时间
fs=20*f;            %Hz       %采样频率
d=1/fs;             %s        %采样间隔
 
 
t=-T/2:d:T/2;       %离散时间t
s1=10*sin(2*pi*10*t+pi/4);    %正弦信号
 
figure(1)
plot(t,s1);
xlabel('时间/s');
ylabel('幅度');

15、符号函数和绘制

  注意，ezplot无法进行分段函数（用function定义）绘制。

lamda=10;
Q=[0,0.5,1,2,10];

syms f;
M1=2*f;
ezplot(M1,[0,2,0,3])

  如果要加颜色和样式，则

ezplot(y1,[0,2,0,3])
% set(findobj('Type','line'),'LineStyle','--','color','r') %注意！会修改这句话前面所有ezplot的绘制样式

参考链接：
[1]Matlab——sym和syms
[2]matlab中的ezplot函数

16、图例等

如果没加figure，那么legend的颜色可能对应不上，出错的代码如下

clc;
clear

lamda=30;
Q=[0,0.5,1,2,10];

syms f;
% figure % 不加这句话，画出来的图的legend颜色不对应
hold on
% y1=1/sqrt((1-1/f^2)/lamda+1);
% ezplot(y1,[0,2,0,3])
% set(findobj('Type','line'),'LineStyle','--','color','r') %会修改这句话前面所有ezplot的绘制样式，所以放最前面
% set(findobj('Type','line'),'LineStyle','--')
M1=1/sqrt(power((1-1/f^2)*Q(1)*f,2)+power((1-1/f^2)/lamda+1,2));
ezplot(M1,[0,2,0,3])
M2=1/sqrt(power((1-1/f^2)*Q(2)*f,2)+power((1-1/f^2)/lamda+1,2));
ezplot(M2,[0,2,0,3])
M3=1/sqrt(power((1-1/f^2)*Q(3)*f,2)+power((1-1/f^2)/lamda+1,2));
ezplot(M3,[0,2,0,3])
M4=1/sqrt(power((1-1/f^2)*Q(4)*f,2)+power((1-1/f^2)/lamda+1,2));
ezplot(M4,[0,2,0,3])
M5=1/sqrt(power((1-1/f^2)*Q(5)*f,2)+power((1-1/f^2)/lamda+1,2));
ezplot(M5,[0,2,0,3])
hold off
title('LLC Gain Cure at lamda=10')
xlabel('f')
ylabel('M(f)')
legend('Q=0','Q=0.5','Q=1','Q=2','Q=10')

17、maltab无法画特别陡的无穷大

  要画的函数如下，根据分母大于0的限制条件可知fs*>0.30151134457776363。

  matlab里画一下如图，可以看到曲线在fs*=0.3233处断了。改大纵轴还是一样。

  代码如下：

clc;
clear

lamda=10;

syms f; % 归一化频率
figure
grid on
hold on
y1=1/sqrt((1-1/f^2)/lamda+1);
ezplot(y1,[0,1,0,10]) % 横轴0-1，纵轴0-10
set(findobj('Type','line'),'LineStyle','--')

18、正弦曲线

figure
grid on
hold on
T=2; % 周期为2
t=0:0.01:1*T; % t从0到2，即画一个周期
w=2*pi/T; 
y=sin(w*t);
plot(t,y)   % 以t为横坐标，则周期为T
plot(w*t,y) % 以wt为横坐标，则周期为2pi
legend('以t为横坐标（T=2）','以wt为横坐标（T=2*pi）')
hold off

19、读取csv后的plot

读取50000x3的csv文件到变量data，然后绘图

plot(data(:,1), data(:,2))

然后报错：tabular/plot 输入参数太多，应该把“()”改为“{}”。

plot(data{:,1}, data{:,2})

20、读取csv文件

第0行开始，第0列开始。
读取整个csv文件：
法一：

M = csvread('data.csv')

法二：读第10行以下，第1列往右

M = csvread('data.csv',10,1)

法三：读第10行到第12行，第1列到第8列之间的数据

M = csvread('data.csv',10,1,[10,1,12,8])

21、fft

fs=100;Ts=1/fs;N=128;   %采样频率和数据点数
n=0:N-1;t=n*Ts;   %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %x为采样得到的信号
y=fft(x,N); % 对采样信号进行N点快速Fourier变换
mag=abs(y); % 求得Fourier变换后的振幅
f=n*fs/N;  % 频率序列（从0到fs）
subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率（即fs/2）之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128 真实幅值');grid on;
subplot(2,2,3),plot(f,[mag(1)/N mag(2:N)*2/N]); %绘出随频率变化的真实振幅，0频率处不需要乘2
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,4)
plot(f(1:N/2),[mag(1)/N mag(2:N/2)*2/N]); %绘出Nyquist频率（即fs/2）之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128 真实幅值');grid on;

参考链接：
[1]MATLAB中fft的使用
[2]Matlab中fft函数的使用与原理
[3]傅里叶变换

22、低通滤波器设计

参考链接：
[1]基于MATLAB的filter的使用，低通、带通和高通滤波器设计
[2]MATLAB Filter Designer实际应用

23、求解任意方程组

用solve，代码如下：

% 单一解
syms x y z
eqn1=x==1;
eqn2=x+y==1;
eqn3=2*x+6*y+3*z==0
[solx,soly,solz]=solve(eqn1,eqn2,eqn3,x,y,z)
solutions=[solx,soly,solz]
% 多个解
syms a u v
[sola,solu,solv]=solve(a*u^2+v^2==0,u-v==1,a^2+6==5*a,a,u,v)
solutions=[sola,solu,solv]

参考链接

24、查询某个值在向量中的位置

idx=find(a==1); % 在向量a中找元素1的位置

25、.m文件中的打印

.m文件中每个语句会以分号结束，如果不以分号结束，就会打印出来。下面代码会打印c和d的值。

a=1;
b=2;
c=3
d=4

24、图注加下标

偶然发现的，但是没法用，因为单一下标用不了。

legend("Gvw_new"); % 可以
legend("Gvw_n") % 不可以

25、化简传函的多项式

num=[1,2];
den=[1,2+3,2*3];
G=tf(num,den); % G=(s+2)/[(s+2)(s+3)];
[num1,den1]=minreal(G.Numerator{:},G.Denominator{:});
Gnew=tf(num1,den1); % Gnew=1/(s+3)

26、取实部、虚部，取模

a=real(complex);
b=imag(complex);
mod=abs(complex);

27、看懂零极点图

  每个极点上的damping指的是：（如果这个极点是某个振荡环节的极点）振荡环节的zeta；frequency指的是振荡环节的wn，overshoot是值这个单一振荡环节的过冲，是根据zeta算出来的。那么问题来了，振荡环节的极点是成对出现的，如果是惯性环节的极点呢？两个惯性环节级联就是振荡环节。
  零点同理，二阶微分环节跟振荡环节对偶。(二阶微分环节是无法画阶跃响应图的，因为是non-casual模型，即分子阶数高于分母）
  那一条红色的曲线是以某个值为半径的圆。

28、关闭注释换行

“主页”-“预设项”-“编辑器/调试器”-“显示”-“右侧文本限制”-取消勾选“显示线条”

29、legend随变量变化为变化

调用num2str()把变量变成字符串，再用strcat()进行字符串拼接。

legend(strcat('Vimin=',num2str(pwr.Vimin),'V/fsmin=',num2str(round(pwm.fsmin/1000)),'kHz/RLmin=',num2str(pwr.RLmin),'ohm'), ...
    strcat('Vinom=',num2str(pwr.Vinom),'V/fr=',num2str(round(pwm.fr/1000)),'kHz/RLmin=',num2str(pwr.RLmin),'ohm'), ...
    strcat('Vimax=',num2str(pwr.Vimax),'V/fsmax=',num2str(round(pwm.fsmax/1000)),'kHz/RLmax=',num2str(pwr.RLmax),'ohm'));

30、读txt文件

调用importdata()函数、fopen()搭配fscanf()函数。

fid=fopen('data.txt','r');
str=fscanf(fid,'%s');
fclose(fid);

31、正弦曲线2

f_sine=50;      % 50Hz正弦波
amp_sine=1;     % 幅值
phase_sine=0; 	% 相位（°）

fs=512e2;       % 采样频率，每隔Δt=1/fs采样一次
N_T=fs/f_sine;  % 采样一个完整周期需要的点数
N=1*N_T;        % 采样1个周期

t_delta=1/fs;
t=0:t_delta:(N-1)*t_delta;
y_sine=amp_sine*sin(2*pi*f_sine*t+pi*phase_sine/180); % 正弦函数为y(t)=Asin(w*t+theta)，获得序列y(0)、y(Δt)、y(2Δt)、y(3Δt)……

32、fft

N=1024;			% FFT点数
fs=100e3;		% 采样频率
f_delta=fs/N;	% 频率分辨率
res_fft=fft(y,N);	% ref_fft为虚数向量

real_amp=abs(res_fft)/(N/2);	% 真正的幅值=模/(N/2)
real_amp(1)=real_amp(1)/2;		% 直流分量=模/N
f=(0:N-1)*f_delta;
figure;
% plot(f,real_amp);	% 横轴为0-fs
plot(f(1:N/2),real_amp(1:N/2));	% 横轴为0-fs/2
title('幅值-频率图');
xlabel('频率');
ylabel('幅值');

real_phase=(1:N);
for i=1:N
	real_phase(i)=phase(res_fft(i));	% 真正的相位（单位为°）
	real_phase(i)=real_phase(i)*180/pi;
end
figure;
% plot(f,real_phase);
plot(f(1:N/2),real_phase(1:N/2));
title('相位-频率图');
xlabel('频率');
ylabel('相位');

33、画对数坐标

x=0:10:1000;
y=0:0.01:1;
semilogx(x,y);	% x轴为对数坐标
semilogy(x,y);	% y轴为对数坐标
loglog(x,y);	% x，y轴都为对数坐标

34、plot的基本要素

xlabel('横轴名称');
ylabel('纵轴名称');
legend('图例1','图例2');
xlim([1 10e3]);	% 设置x轴范围
set(gca,'XTick',[1:2:10e3]); % 设置x轴范围，以及x轴间隔。y轴类似
axis([1 10e3 2 20e3]);	% 同时控制x轴y轴范围
title('图名');

35、脚本文件无法复制粘贴

原因：初始化的Java堆内存不够，在设置中增加内存即可。
设置方法：预设（Preferences）-> MATLAB -> 常规（General） -> Java堆内存（Java Heap Memory），利用滑块或数字增大Java堆内存即可。

36、导出伯德图的数据

% 构造频率采样点，用于导出伯德图数据
% 下面构造的w_matrix为
% 10e-4 11e-4 ... 99e-4	（Hz）
% 100e-4 110e-4 ... 990e-4
% ...
% 1e7 ...	99e7
w_matrix=zeros(11,90);
w=zeros(1,990);
for i=1:size(w_matrix,1)
    for j=1:size(w_matrix,2)
        w_matrix(i,j)=2*pi*(10^(i-1)*10*1e-4)+(j-1)*10^(i-1)*1e-4;
        w((i-1)*size(w_matrix,2)+1:(i-1)*size(w_matrix,2)+size(w_matrix,2))=w_matrix(i,:);
    end
end

[mag,phase]=bode(Gs,w);	% 导出伯德图数据，这里调用bode不会产生图像
mag_new=reshape(mag(1,1,:),1,[]);	% mag和phase都是3维数据，要reshape成一维，注意，这里的mag是A，如果要用db表示，需要20*log10(mag)
phase_new=reshape(phase(1,1,:),1,[]);

37、对伯德图导出来的数据进行phase warp

function phase_new=phase_warp(phase)
phase_new=phase;
% 排除维度影响
if(size(phase_new,1)==1)
    len=size(phase_new,2);
else
    len=size(phase_new,1);
end
% warp
for i=1:len
    if(phase_new(i)>180)
        while(phase_new(i)>180)
            phase_new(i)=phase_new(i)-360;
        end
    elseif(phase_new(i)<-180)
        while(phase_new(i)<-180)
            phase_new(i)=phase_new(i)+360;
        end
    end
end
end

38、读取excel

sheet=xlsread('data.xlsx');

39、伯德图设置

opts=bodeoptions;
opts.xlim={[1e-2 10e6]};
opts.freqUnits='Hz';
opts.PhaseWrapping='on';

figure;
bode(Gs,opts);

40、传递函数构造

以一阶惯性环节为例

fc=1e3; % 截止频率为1kHz
num=1；
den=[1/(2*pi*fc) 1];
figure;
G=tf(num,den);
title('一阶惯性环节（fc=1kHz）');
grid on

41、求某个频率下的正弦波延时

已知正弦波频率为f，单位Hz，该频率下的相移为phase，则

tao=(phase*pi/180)/(2*pi*f);

42、两个传递函数级联

G=series(G1,G2);

43、（传递函数）两个多项式相乘

num3=conv([1/sys_c3.wz1,1],[1/sys_c3.wz2,1]);
den3=conv([1,0],[1/sys_c3.wp,1]);
sys_c3.Gc=tf(num3,den3);

44、提取高维数组的一个维度

用reshape

array=round(10*randn(2,3,4)); % 一个2x3x4维的三维数组
data=reshape(array(:,1,1),1,[]); % 取坐标为(X,1,1)的数，总共2个

45、plot线条/标记加粗

plot(x, y, 'b-*', 'Linewidth', 2, 'MarkerSize', 8); % Linewidth是线条加粗、MarkerSize是标记加粗