matlab仿真介绍,基于MATLAB仿真频率细化过程介绍

1 频率细化过程介绍

频率细化是在信号处理和模态分析中广泛应用的一种技术,它能够提高频率的分辨率,将选定的频率域上的特性曲线放大,从而使系统的频率特性能更清楚地显示出来,如图1所示。

70b873ca61b0a9cb6e7e219985742159.gif

设系统的采样频率为fs,采样点数为NO,则频率分辨率为:

Δf=fs/NO

从上式可以看出,要进行频率细化,即提高频率分辨率,使Δf变小,有两种方法:增加采样点数NO和降低采样频率fs,这里只介绍降低采样频率的方法。

这种方法主要是基于移频原理,如图2所示。

633f6bebc18638547e112b079369d311.gif

设想要移频部分的频率为fp,其角频率为wp=2πfp,令f (t)=exp(-jwpt)。时域信号x(t)与f(t)卷积后,则在频域上,该信号的fp频率就移到了原点处。

信号频率移到低频后,经过低通滤波,就可以用低的采样频率进行采样,从而达到提高频率分辨率、频率细化的目的。具体过程如图3所示。

a5d1cc4eab43faf2655ca3429d3ae960.gif

2 移频法频率细化过程中几个问题的分析

(1) 非细化处理的系统采样频率为fs,采样点数为NO。对于细化过程,设频率细化倍数为Nr,信号经过抗混叠滤波器后进行A/D采样,采样频率应仍为原来的 fs,保持不变,采样点数则为Nr*NO,这样就保证了细化与非细化处理的基本频带范围保持不变,并且可以细化这一频带中的任何一段。

(2)要细化的频率范围为fl~fu,移频后,低频点fl移到原点,则高频点变为(fu-fl)。这时数字低通滤波器的截止频率应大于(fu-fl),并小于低频重采样频率fs/Nr的一半。即截止频率的范围为:

(fu-fl)<fc<fs/2Nr

且可以得到最大细化倍数N与细化频率范围之间的关系为:

Nr<fs/2(fu-fl)

这为设置细化倍数范围提供了依据。

ce4dc92869b76159db05e6a3b2713872.gif

(3)数字低通滤波器的通带必须平,通带内波动要小,这样原信号的频率特性细化后在幅值上才不会改变;同时,最好使滤波器的带外衰减>-70dB,且-70dB处的频率<fs/2Nr,这样就能保证低频重新采样时抗混叠的效果最好,细化的效果也最好。

(4) 细化与不细化过程占用时间的比较:由于采样点数NO保持不变,因此细化处理的FFT时间与不细化的FFT时间一样,都是NOLog(NO)/2;细化过程 要进行Nr*NO点的高频采样和NO点的低频采样,而不细化过程只进行NO点的高频采样,所以在采样时间上,细化过程要稍长一点。但它与NO点的FFT变 换时间比起来可以不计,因此,细化处理与不细化处理在时间上差不多。

e42be0dd2a688643b2ba18847f700e14.gif

(5) 移频法频率细化与增加采样点数频率细化的比较:移频法频率细化只进行NO个点的FFT变换,和一些数组、矩阵的运算,它所花的时间约 为:NOLog(NO)/2;而采用增加采样点数频率细化要进行Nr*NO个点的FFT变换,它所花的时间约为:(Nr*NO)Log(Nr*No) /2,因此采用移频法频率细化的时间要短得多,这是它的优点,也是工程应用中多采用这种方法的原因,我们也只对它进行仿真。

由于移频使fl前面的频段移到频域的负轴上,而低通滤波又滤掉了fu后面的高频部分,因此,这种方法只能进行一段频率的细化,不能进行全频段的细化,这是移频法频率细化的缺点。要进行全频段细化,可以采用增加采样点数的方法。

bc6dc76d77fb660f651127755fc3f26d.gif

3 用MATLAB仿真频率细化过程

用MATLAB程序仿真图3这个过程,主要实现A/D采样、移频、低通滤波、低频重新采样、FFT变换等,同时注意上面几个问题的分析。

总 信号由302Hz、304Hz、306Hz、308Hz四个不同频率、不同幅值的正余弦信号合成。采样点为512,采样频率为5120Hz,则频率分辨率 为10Hz,这在频域内分辨不出这四个信号。仿真软件在300~320Hz范围内细化10倍,则这时频率分辨率为1Hz,就可以逐渐看清楚这四个信号。

25d1b796da6175a18f74d87234be524e.gif

x(t)=sin(2*302)+2cos(2π*304)+3sin(2π*306)+4sin(2π*308)

细化处理需要5120个采样点,仿真时A/D采样用对连续信号x(t) 进行离散化来代替,离散化的点值就是采样值,它们组成一个数组(矩阵),这是一个点序列。

采样完成后,对采样点进行移频仿真(将300Hz移频到原点),由以下步骤实现:

令:F1=300W1=2π*F1

f(t)=exp(-jw1*t)

t取离散化时刻,Δt=1/fsTk=k*Δt

则f(t)离散化可变为:

f(k)=exp(-jw1*Tk)

=exp[(-j2π*f1*k)/fs]

这也是一个点序列。

MATLAB将上面两个点序列进行卷积(矩阵相乘),得到一个新的点序列,就实现了移频过程,即使F1=300Hz 在频域上移动了0,其它频率特性依此前移。

MATLAB语言有很方便的滤波器设计和数字滤波功能,可用它设计一个带内波动<0.1dB,带外衰减>-70dB,截止频率为20Hz的低通滤波器,并将上面的点序列进行数字滤波,得到一个只在0~20Hz频段上有特性曲线、其它频段被滤掉的信号的点序列。

对上面的点序列每隔10点(细化倍数)进行抽取,实际上就是对经移频、滤波后的信号进行重新采样。采样频率为S120Hz/10,即采样频率降低了10倍,抽取得到512个点。

0dcf63a63ae2f2fb444b4950c4ae4112.gif

对这512个点进行FFT变换,就可以得到300Hz~320Hz频段上已细化10倍的频率特性曲线了,结果如图4、图8所示。

用仿真程序分别进行2、4、8、10倍的频率细化,从运行结果可以看出,细化倍数越大,频率分辨率就越高,频域波形越看得清楚,这与理论是相符合的。

MATLAB中使用FFT做频谱分析时频率分辨率问题-频率分辨率.rar MATLAB中使用FFT做频谱分析时频率分辨率问题 最近做FFT时,使用的采样频率和信号长度的取舍一直没有搞清楚,后来在论坛上发了一个贴子《总结一下使用FFT和维纳-辛钦定理求解PSD问题》(讨论见https://www.ilovematlab.cn/thread-27150-1-1.html,特别感谢会员songzy41,他的问题给了我很大启示),跟帖中给了我不少启示,并且让我对“频率分辨率”这个概念有了更深入的理解。再次一并感谢论坛的高手们。 频率分辨率,顾名思义,就是将信号中两个靠的很近的频谱分开的能力。 信号x长度为Ts,通过傅氏变换后得到X,其频率分辨率为Δf=1/T(Hz),若经过采样后,假设采样频率为fs=1/Ts,而进行频谱分析时要将这个无穷长的序列使用窗函数截断处理,假设使用矩形窗,我们知道,矩形窗的频谱为sinc函数,主瓣宽度可以定义为2*pi/M,M为窗宽,那么,时域相乘相当于频域卷积,频域内,这一窗函数能够分辨出的最近频率肯定不可能小于2*pi/M了,也就是如果数据长度不能满足2*pi/M<|w2-w1|(w2,w1为两个靠的很近的频率),那么在频谱分析时,频谱上将不能分辨出这两个谱,由于w2-w1=2*pi/fs=2*pi*Δf/fs也就是2*pi/M<2*piΔf/fs,得到Δf的限制为fs/M,这就是窗函数宽度的最小选择,就是说,根据Shannon采样定理确定了采样频率后,要根据靠的最近的谱峰来确定最小的采样长度,这样,所作出来的频谱才能分辨出那两个谱峰,也就是拥有了相应的频率分辨率。 几个例子: 考虑双正弦信号:x = sin sin;根据Shannon采样定理,采样频率要大于截止频率的两倍,这里选采样频率为80,那么,我们可以看到,Δf为0.2Hz,那么,最小的数据长度为0.2/80=400,但是对正弦信号的频谱分析经验告诉我们,在截断时截断时的数据要包含整周期,并且后面不宜补零以避免频谱泄露(这一点见胡广书《数字信号处理导论》,清华大学出版社),那么,我们要选择至少980个点,才能保含到一个整周期,另外,FFT的经验告诉我们作分析时最好选择2的整数次幂,我们选择靠的最近的1024点。分析结束。 [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; figure; X = fftshift); plot*Fs/N,abs*2/N); grid on; axis; 这是按照我们的分析进行的编程和图形 zheng.jpg 可以看出这两个谱峰很好的被分辨开来,9.8Hz不在谱线上,所以幅值不为1,以下是一些对比: [CODE] Fs = 80; n = 0:1/Fs:1023*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:979*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; n = 0:1/Fs:399*1/Fs; x = sin sin; N = length; X = fftshift); figure; subplot plot*Fs/N,abs*2/N); grid on; axis; title; Fs = 20; n = 0:1/Fs:1024*1/Fs; x = sin sin; N = length; X = fftshift); subplot plot*Fs/N,abs*2/N); grid on; axis; title; 结果如下: 1024.jpg 400.jpg 这是我在做FFT以及论坛中的问题时所得到的一点启发,不当之处还请大家指正。OO~ 频率分辨率.rar 为了方便大家,我将doc版报告和m文件一起上传,和帖子内容一样。OO~
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值