查看: 6495|回复: 2|关注: 0
求教：Matlab在pid的微分滤波器模块中，系数N作用？
<h1 style="color:# 麦片财富积分
新手, 积分 5, 距离下一级还需 45 积分
求教： 下面的pid中的微分滤波器有两个作用： 1，差分& &2，滤波
请教各位， Filter coefficient N 是怎么影响滤波效果的？
该模块的simulink框图和传递函数如下图：
<h1 style="color:# 麦片财富积分
这样很像个比例系数。
<h1 style="color:# 麦片财富积分
回复 2# tangajun 的帖子
是的。看着像。
这个传递函数是高通的
站长推荐 /2
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区查看: 42038|回复: 37|关注: 0
MATLAB信号分析中对频率谱、相位谱、功率谱及相关性的处理
<h1 style="color:# 麦片财富积分
新手, 积分 8, 距离下一级还需 42 积分
关注者： 8
& && &声明：这也是本人搜集的，不是出于本人之手。不好请见谅，或是引了用某位大神的帖子，请原谅！
& && && &&&MATLAB处理信号得到频谱、相谱、功率谱
第一：频谱
一.调用方法
X=FFT(x)；
X=FFT(x，N)；
x=IFFT(X);
x=IFFT(X,N)
用MATLAB进行谱分析时注意：
（1）函数FFT返回值的数据结构具有对称性。
xn=[4 3 2 6 7 8 9 0];
Xk=fft(xn)
39.0000& && && && &-10.7782 + 6.2929i& && && &0 - 5.0000i& & 4.7782 - 7.7071i& & 5.0000& && && && &&&4.7782 + 7.7071i& && && &0 + 5.0000i -10.7782 - 6.2929i
Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。
（2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。
二.FFT应用举例
例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
fs=100;N=128;& & %采样频率和数据点数
n=0:N-1;t=n/& & %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N);& &&&%对信号进行快速Fourier变换
mag=abs(y);& && &%求得Fourier变换后的振幅
f=n*fs/N;& &&&%频率序列
subplot(2,2,1),plot(f,mag);& & %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');
%对信号采样数据为1024点的处理
fs=100;N=1024;n=0:N-1;t=n/
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N);& & %对信号进行快速Fourier变换
mag=abs(y);& & %求取Fourier变换的振幅
subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');
运行结果：
& & & & & & & & & & & & & & & &
& && &&&fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。
例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：
（1）数据个数N=32，FFT所用的采样点数NFFT=32；
（2）N=32，NFFT=128；
（3）N=136，NFFT=128；
（4）N=136，NFFT=512。
fs=100; %采样频率
Ndata=32; %数据长度
N=32; %FFT的数据长度
n=0:Ndata-1;t=n/& & %数据对应的时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);& & %时间域信号
y=fft(x,N);& & %信号的Fourier变换
mag=abs(y);& &&&%求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');
Ndata=32;& & %数据个数
N=128;& && &%FFT采用的数据长度
n=0:Ndata-1;t=n/& & %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');
Ndata=136;& & %数据个数
N=128;& && &%FFT采用的数据个数
n=0:Ndata-1;t=n/ %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;& & %真实频率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');
Ndata=136;& &&&%数据个数
N=512;& &&&%FFT所用的数据个数
n=0:Ndata-1;t=n/ %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;& & %真实频率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅
xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');
& & & & & & & & & & & & & & & &
（1）当数据个数和FFT采用的数据个数均为32时，频率分辨率较低，但没有由于添零而导致的其他频率成分。
（2）由于在时间域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。
（3）FFT程序将数据截断，这时分辨率较高。
（4）也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT振幅谱也基本不受影响。
& && &对信号进行频谱分析时，数据样本应有足够的长度，一般FFT程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。
例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)
& & & & & & & & & & & & & & & &
（1）数据点过少，几乎无法看出有关信号频谱的详细信息；
（2）中间的图是将x(n)补90个零，幅度频谱的数据相当密，称为高密度频谱图。但从图中很难看出信号的频谱成分。
（3）信号的有效数据很长，可以清楚地看出信号的频率成分，一个是0.24Hz，一个是0.26Hz，称为高分辨率频谱。
& && && &可见，采样数据过少，运用FFT变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数，谱的密度增高了，但仍不能分辨其中的频率成分，即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。
第二： 相谱
（相位谱和频率普是回事儿，想着把频谱中的幅值部分换成相角就可以了）
&&由于没有找到具体的理论，就举几个例子说明一下。
& & 比如要求y=sin(2*pi*60*t) 的相位谱，
程序如下：
fs=200;N=1024;n=0:N-1;t=n/y=sin(2*pi*60*t);
Y=fft(y,N);
A=abs(Y);f=n*fs/N;
ph=2*angle(Y(1:N/2));
ph=ph*180/
plot(f(1:N/2),ph(1:N/2));
xlabel('频率/hz'),ylabel('相角'),title('相位谱');
期中的 ph=2*angle(Y(1:N/2));ph=ph*180/是利用angle函数求出每个点的角度，并由弧度转化成角度！
angle函数解释：
Phase angle
P = angle(Z)
Description
P = angle(Z) returns the phase angles, in radians, for each element of complex array Z. The angles lie between ±π.
For complex Z, the magnitude R and phase angle theta are given by
R = abs(Z)
theta = angle(Z)
and the statement
Z = R.*exp(i*theta)
converts back to the original complex Z.
Z = [ 1 - 1i& &2 + 1i& &3 - 1i& &4 + 1i
& && &1 + 2i& &2 - 2i& &3 + 2i& &4 - 2i
& && &1 - 3i& &2 + 3i& &3 - 3i& &4 + 3i
P = angle(Z)
& &-0.7854& & 0.4636& &-0.3218& & 0.2450
& & 1.1071& &-0.7854& & 0.5880& &-0.4636
& &-1.2490& & 0.9828& &-0.7854& & 0.6435
& & 1.3258& &-1.1071& & 0.9273& &-0.7854
Algorithms
The angle function can be expressed as angle(z) = imag(log(z)) = atan2(imag(z),real(z)).
第三：功率谱
matlab实现经典功率谱估计
fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数
matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch法估计，即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。
1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码示例：
Fs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
window=boxcar(length(xn)); %矩形窗
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
plot(f,10*log10(Pxx));
2、间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
Matlab代码示例：
Fs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.1、Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
Matlab代码示例：
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(length(n)); %矩形窗
noverlap=0; %数据无重叠
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
plot(k,plot_Pxx);
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
3.2、Welch法
Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Matlab代码示例：
n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
nfft=1024;
window=boxcar(100); %矩形窗
window1=hamming(100); %海明窗
window2=blackman(100); %blackman窗
noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率
[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
plot(f,plot_Pxx);
plot(f,plot_Pxx1);
plot(f,plot_Pxx2);
第四： 相关性分析
1. 首先说说自相关和互相关的概念。& & 这个是信号分析里的概念，他们分别表示的是两个时间序列之间和同一个时间序列在任意两个不同时刻的取值之间的相关程度，即互相关函数是描述随机信号x(t),y(t)在任意两个不同时刻t1，t2的取值之间的相关程度，自相关函数是描述随机信号x(t)在任意两个不同时刻t1，t2的取值之间的相关程度。& & 自相关函数是描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个判断指标，把两测点之间信号的互谱与各自的自谱联系了起来。它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生的误差非常有效.& & 事实上，在图象处理中，自相关和互相关函数的定义如下：设原函数是f(t)，则自相关函数定义为R(u)=f(t)*f(-t)，其中*表示卷积；设两个函数分别是f(t)和g(t)，则互相关函数定义为R(u)=f(t)*g(-t)，它反映的是两个函数在不同的相对位置上互相匹配的程度。那么，如何在matlab中实现这两个相关并用图像显示出来呢？dt=.1;
t=[0:dt:100];
[a,b]=xcorr(x,'unbiased');
plot(b*dt,a)
上面代码是求自相关函数并作图，对于互相关函数，稍微修改一下就可以了，即把[a,b]=xcorr(x,'unbiased');改为[a,b]=xcorr(x,y,'unbiased');便可。
2. 实现过程：
& && &在Matalb中，求解xcorr的过程事实上是利用Fourier变换中的卷积定理进行的，即R(u)=ifft(fft(f)×fft(g))，其中×表示乘法，注：此公式仅表示形式计算，并非实际计算所用的公式。当然也可以直接采用卷积进行计算，但是结果会与xcorr的不同。事实上，两者既然有定理保证，那么结果一定是相同的，只是没有用对公式而已。下面是检验两者结果相同的代码：
t=[0:dt:100];
x=3*sin(t);
y=cos(3*t);
subplot(3,1,1);
plot(t,x);
subplot(3,1,2);
plot(t,y);
[a,b]=xcorr(x,y);
subplot(3,1,3);
plot(b*dt,a);
yy=cos(3*fliplr(t)); % or use: yy=fliplr(y);
z=conv(x,yy);
subplot(3,1,3);
plot(b*dt,z,'r');
即在xcorr中不使用scaling。3. 其他相关问题：
（1）相关程度与相关函数的取值有什么联系？& & 相关系数只是一个比率，不是等单位量度，无什么单位名称，也不是相关的百分数，一般取小数点后两位来表示。相关系数的正负号只表示相关的方向，绝对值表示相关的程度。因为不是等单位的度量，因而不能说相关系数0.7是0.35两倍，只能说相关系数为0.7的二列变量相关程度比相关系数为0.35的二列变量相关程度更为密切和更高。也不能说相关系数从0.70到0.80与相关系数从0.30到0.40增加的程度一样大。
对于相关系数的大小所表示的意义目前在统计学界尚不一致，但通常按下是这样认为的：
相关系数& && &相关程度
0.00-±0.30& & 微相关
±0.30-±0.50&&实相关
±0.50-±0.80&&显著相关
±0.80-±1.00&&高度相关（2）matlab计算自相关函数autocorr和xcorr有什么不一样的？& & 分别用这两个函数对同一个序列计算，为什么结果不太一样？因为xcorr是没有将均值减掉做的相关，autocorr则是减掉了均值的。而且，用离散信号做自相关时，信号截取长度（采样点N）不一样，自相关函数就不一样。
（3）xcorr是计算互相关函数，带有一个option的参数:
a=xcorr(x,y,'option')
option=baised时，是计算互相关函数的有偏估计；
option=unbaised时，是计算互相关函数的无偏估计；
option=coeff时，是计算归一化的互相关函数，即为互相关系数，在-1至1之间；
option=none，是缺省的情况。
所以想要计算互相关系数，可用'coeff'参数。
*************************************************************************
用这个xcorr函数作离散互相关运算时要注意，当x, y是不等长向量时，短的向量会自动填0与长的对齐，运算结果是行向量还是列向量就与x一样。
互相关运算计算的是x,y两组随机数据的相关程度，使用参数coeff时，结果就是互相关系数，在-1至1之间，否则结果不一定在这范围，有可能很大也有可能很小，这视乎x, y数据的大小，所以一般要计算两组数据的相关程度，一般选择coeff参数，对结果进行归一化。
所谓归一化简单理解就是将数据系列缩放到-1到1范围，正式的就是一种简化计算的方式，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为纯量。变换式为X=(X实测--Xmin)/(Xmax-Xmin)。
一般来说选择归一化进行互相关运算后，得到结果绝对值越大，两组数据相关程度就越高。
<h1 style="color:# 麦片财富积分
好贴，通过几种方法的仿真，更加清楚了解如何分析信号的功率谱，而且更加了解了它们的意义。
<h1 style="color:# 麦片财富积分
关注者： 4
终于找到这样的好帖子了，虽然是收集别人的，但是还是很经典的，赞一下！:lol
<h1 style="color:# 麦片财富积分
谢谢楼主的搜集，把几种方法集合到一起，最近正在学这个，看了感觉清楚了很多
<h1 style="color:# 麦片财富积分
十分感谢楼主，学习了
<h1 style="color:# 麦片财富积分
<h1 style="color:# 麦片财富积分
学习了，多谢楼主
<h1 style="color:# 麦片财富积分
好贴，赞一个
<h1 style="color:# 麦片财富积分
好贴，解决我的问题啊。
<h1 style="color:# 麦片财富积分
确实很有用，谢谢楼主
站长推荐 /2
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区查看: 2961|回复: 5|关注: 0
simulink 中坐标变换的系数为什么是2/3？？
<h1 style="color:# 麦片财富积分
新手, 积分 38, 距离下一级还需 12 积分
关注者： 2
同题，为什么simulink中从三相到两相的静止坐标变换矩阵系数是2/3，而不是sqrt(2/3)。而且变换矩阵的最后一行为什么是1/2？？？？
<h1 style="color:# 麦片财富积分
关注者： 2
在线等，跪求大侠指教！！
<h1 style="color:# 麦片财富积分
关注者： 2
没人啊？求指教啊！！！
另外，再提个问题，就是采用simulink中的坐标变换方法，是否满足功率不变？
<h1 style="color:# 麦片财富积分
这是电机变换折算过来的，三相折算两相，还有运动折算静止，你看一下电机的坐标变换，
<h1 style="color:# 麦片财富积分
关注者： 1
|此回复为最佳答案
如果按照等功率变换，则变换系数为sqrt（2、3）,若为等幅值变换，则系数为2、3.
<h1 style="color:# 麦片财富积分
关注者： 2
highspeed1989 发表于
如果按照等功率变换，则变换系数为sqrt（2、3）,若为等幅值变换，则系数为2、3. ...
假设定转子最大互感是Lms，那么在两相同步旋转坐标系下，互感Lm=1.5Lms(此结果为等功率变换下，推导结果)。那么在估计磁链时，那么磁链应大概等于Lm*但仿真结果却发现磁链等于Lms*ism。
这是为什么呢？
站长推荐 /2
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区