声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

声振论坛 展示 振动理论 信号处理 查看内容

运用 FFT 求取矩形脉冲的频谱

2015-11-3 08:00| 发布者: aspen| 查看: 1116| 评论: 2|原作者: happy|来自: 声振论坛

摘要: % 本程序采用 FFT 计算连续时间 Fourier 变换。输出幅频谱数据对 (f,AW) 。 % 输入量 (wt,t) 为已经窗口化了的时间函数 wt(t) ,它们分别是长度为 N 的向量。 % 对于 " 非平凡 " 取值时段有限的情况,应使该时 ...
运用 FFT 求取矩形脉冲的频谱
[cftbyfft.m]

function [AW,f]=cftbyfft(wt,t,flag)
%cftbyfft.m
% 本程序采用 FFT 计算连续时间 Fourier 变换。输出幅频谱数据对 (f,AW) 。
% 输入量 (wt,t) 为已经窗口化了的时间函数 wt(t) ,它们分别是长度为 N 的向量。
% 对于 " 非平凡 " 取值时段有限的情况,应使该时段与窗口长度相比足够小。以
% 提高频率分辨率。
% 对于 " 非平凡 " 取值时段无限的情况,窗口长度的选取应使窗口外的函数值小
% 到可忽略,以提高近似精度。
% 输入宗量 flag 控制输出 CFT 的频率范围。
% flag 取非 0 时(缺省使用),频率范围在 [0,fs );
% flag 取 0 时,频率范围在 [-fs/2,fs/2) 。
if nargin==2;flag=1; end
N=length(t); % 采样点数,应为 2 的幂次,以求快速。
T=t(length(t))-t(1); % 窗口长度
dt=T/N; % 时间分辨率。
W0=fft(wt); % 施行 FFT 变换 <16>
W=dt*W0; % 算得 [0,fs) 上的 N 点 CFT 值
df=1/T; % 频率分辨率
n=0:1:(N-1);
% 把以上计算结果改写到 [-fs/2,fs/2] 范围
if flag==0
n=-N/2:(N/2-1);
W=fftshift(W); % 产生满足式( 5.13.3.1-6 )的频谱
end

f=n*df; % 频率分度向量
AW=abs(W); % 福频谱数据向量
if nargout==0
plot(f,AW);grid,xlabel( ' 频率 f' );ylabel( '|w(f)|' )
end

运行以下指令,绘制时域波形和幅频谱
M=5; % 做 2 的幂次用。本例把 M 设得较小,是为了观察混迭。 <1>
tend=1; % 波形取非零值的时间长度。
T=10; % 窗口化长度应足够大,以减小窗口化引起的泄露“旁瓣”效应。 <3>
N=2^M; % 采样点数,取 2 的幂是为使 FFT 运算较快。
dt=T/N; % 以上 T 、 N 的取值应使 N/T=fs 采样频率大于两倍时间波形带宽,以克服
% 采样引起的频谱混迭。
% 在本例中,据理论分析知 W(f=7.5)=Sa(7.5*pi)=1/(7.5*pi)<5% 。
% 因此,可近似认为本例时间信号带宽为 7.5Hz 。
n=0:N-1; % 采样序列
t=n*dt; % 采样点时间序列
w=zeros(size(t,2),1);
Tow=find((tend-t)>0); % 产生非零波形时段的相应序列
w(Tow,1)=ones(length(Tow),1); % 在窗口时段内定义的完整波形
plot(t,w,'b','LineWidth',2.5),title('Time Waveform');xlabel('t --- >') 

[AW,f]=cftbyfft(w,t,0); 
ff=f+eps; % 为避免下面指令出现 0/0 而采取的措施
AWW=abs(sin(pi*ff)./(pi*ff));
plot(f,AW,'b-',ff,AWW,'r:')
title('Aliasing caused by undersampling')
xlabel('f --- >');ylabel('|W(f)|'),legend('by FFT','Theoretical') 

如果采样频率太低会出现混迭现象
发表评论

最新评论

引用 happy 2006-3-25 08:15
不知道pwm脉冲波什么样,一个例子

运用 FFT 求取矩形脉冲的频谱

[cftbyfft.m]
  1. function [AW,f]=cftbyfft(wt,t,flag)
  2. %cftbyfft.m
  3. % 本程序采用 FFT 计算连续时间 Fourier 变换。输出幅频谱数据对 (f,AW) 。
  4. % 输入量 (wt,t) 为已经窗口化了的时间函数 wt(t) ,它们分别是长度为 N 的向量。
  5. % 对于 " 非平凡 " 取值时段有限的情况,应使该时段与窗口长度相比足够小。以
  6. % 提高频率分辨率。
  7. % 对于 " 非平凡 " 取值时段无限的情况,窗口长度的选取应使窗口外的函数值小
  8. % 到可忽略,以提高近似精度。
  9. % 输入宗量 flag 控制输出 CFT 的频率范围。
  10. % flag 取非 0 时(缺省使用),频率范围在 [0,fs );
  11. % flag 取 0 时,频率范围在 [-fs/2,fs/2) 。
  12. if nargin==2;flag=1; end
  13. N=length(t); % 采样点数,应为 2 的幂次,以求快速。
  14. T=t(length(t))-t(1); % 窗口长度
  15. dt=T/N; % 时间分辨率。
  16. W0=fft(wt); % 施行 FFT 变换 <16>
  17. W=dt*W0; % 算得 [0,fs) 上的 N 点 CFT 值
  18. df=1/T; % 频率分辨率
  19. n=0:1:(N-1);
  20. % 把以上计算结果改写到 [-fs/2,fs/2] 范围
  21. if flag==0
  22. n=-N/2:(N/2-1);
  23. W=fftshift(W); % 产生满足式( 5.13.3.1-6 )的频谱
  24. end

  25. f=n*df; % 频率分度向量
  26. AW=abs(W); % 福频谱数据向量
  27. if nargout==0
  28. plot(f,AW);grid,xlabel( ' 频率 f' );ylabel( '|w(f)|' )
  29. end

  30. 运行以下指令,绘制时域波形和幅频谱
  31. M=5; % 做 2 的幂次用。本例把 M 设得较小,是为了观察混迭。 <1>
  32. tend=1; % 波形取非零值的时间长度。
  33. T=10; % 窗口化长度应足够大,以减小窗口化引起的泄露“旁瓣”效应。 <3>
  34. N=2^M; % 采样点数,取 2 的幂是为使 FFT 运算较快。
  35. dt=T/N; % 以上 T 、 N 的取值应使 N/T=fs 采样频率大于两倍时间波形带宽,以克服
  36. % 采样引起的频谱混迭。
  37. % 在本例中,据理论分析知 W(f=7.5)=Sa(7.5*pi)=1/(7.5*pi)<5% 。
  38. % 因此,可近似认为本例时间信号带宽为 7.5Hz 。
  39. n=0:N-1; % 采样序列
  40. t=n*dt; % 采样点时间序列
  41. w=zeros(size(t,2),1);
  42. Tow=find((tend-t)>0); % 产生非零波形时段的相应序列
  43. w(Tow,1)=ones(length(Tow),1); % 在窗口时段内定义的完整波形
  44. plot(t,w,'b','LineWidth',2.5),title('Time Waveform');xlabel('t --- >')

  45. [AW,f]=cftbyfft(w,t,0);
  46. ff=f+eps; % 为避免下面指令出现 0/0 而采取的措施
  47. AWW=abs(sin(pi*ff)./(pi*ff));
  48. plot(f,AW,'b-',ff,AWW,'r:')
  49. title('Aliasing caused by undersampling')
  50. xlabel('f --- >');ylabel('|W(f)|'),legend('by FFT','Theoretical')
复制代码



如果采样频率太低会出现混迭现象
引用 cycr1234 2006-3-25 11:47
不好意思错了~~
叫付利叶变换叫习惯了
是对无限周期函数求付利叶极数

查看全部评论(2)

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 09:53 , Processed in 0.044752 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部