基于matlab的矩形波表⽰,对MATLAB中傅⽴叶函数的解读-以
矩形波的傅⾥叶变换为例...
题⽬:对⼀个矩形波进⾏傅⽴叶变换,并观察其频谱与理想⽅波函数的理论频谱是否相同
代码如下:
Fs=100; % Sampling frequency
T=1/Fs; % Sample time
L=100; % Length of signal
t=(0:L-1)*T; % Time vector
%绘制占空⽐为40%的矩形波
y=zeros(1,L); y(1:0.4*L)=1;
subplot(121); plot(t,y);
title('占空⽐为40%的矩形波');
xlabel('Time'); ylabel('Amplitude');
%对矩形波进⾏傅⽴叶变换
NFFT=1024;
Y=fft(y,NFFT);
f=Fs*linspace(0,1,NFFT);
subplot(122); plot(f,abs(Y(1:NFFT))) ;
title('占空⽐为40%的矩形波对应的频谱曲线');
xlabel('Frequency (Hz)'); ylabel('|Y(f)|');
代码运⾏结果:
图1
编程解析:
对于周期性⽅波⽽⾔,其各个周期信号的频谱会产⽣重叠,不便于观察,因此,在程序中未使⽤square(
)函数,⽽是通过使⽤zeros( )函数,然后对其进⾏部分置1来构造⼀个矩形波。
对矩形波进⾏傅⽴叶变换部分,其中
Y=fft(y,NFFT)
表⽰对时域波形进⾏NFFT点的离散傅⽴叶变化,NFFT数值越⼤得到的频域图形其曲线越平滑,这⾥的NFFT最好取为2的整数次幂,因为当NFFT为2的整数次幂的情况下进⾏离散傅⽴叶变换需要的次数就越少,效率就越⾼。⽽⼀般情况下,NFFT应该取⼤于信号长度L且距L最近的⼀个2的整数次幂,在MATLAB中⽤NFFT=2^nextpow2(L)实现。
当进⾏离散傅⽴叶变换的点数⼤于信号的长度L时,会对L后⾯的部分补零,当NFFT的数值⼩于L时,会取L的前NFFT个点,⽽把后⾯的点舍弃。在进⾏傅⽴叶变换时,倘若不使⽤NFFT,⽽直接写成fft(y)的形式,此时是默认对时域信号y进⾏L点的离散傅⽴叶变换。
⽽ f=Fs*linspace(0,1,NFFT) 表⽰在0到Fs之间取NFFT个点,作为频谱曲线的横轴。在MATLAB
help中给出的例⼦为:
NFFT=2^nextpow2(L); %NFFT为距离L最近的2的整数次幂的⼀个数
Y=fft(y,NFFT)/L; %矫正频谱信号的幅值
f=Fs/2*linspace(0,1,NFFT/2+1); %此处作的是单边频谱
plot(f,2*abs(Y(1:NFFT/2+1))); %对于单边功率谱,其幅值应该乘以2
以上例⼦是仅仅作出时域信号的单边频谱,因此定义频谱的横轴时只取了NFFT的⼀半
即linspace(0,1,NFFT/2+1) 。此时得到的频谱横轴还仅仅是归⼀化的数值,MATLAB中规定对于单边频谱需要将其归⼀化的横轴乘以采样频率Fs的⼀半,即Fs/2。对于双边频谱,应将其归⼀化横轴乘以Fs,此时的得到的信号的时域频域才是完美对应的。
仿照MATLAB
help 中给出的例⼦,画出占空⽐为40%的矩形波的时域频谱图像如下所⽰:
图2
很容易发现此时信号的频谱图只有Fs的⼀半,但是频域信号的第⼀个零点仍然为1/T,符合理论情况。同时可以看出频谱信号的幅值都变为原来的⼆倍。
倘若将
f=Fs/2*linspace(0,1,NFFT/2+1)中的Fs/2去除,或者乘以其他的数值,这样的话信号的频域时域信号便不能和理论情况对应。下⾯是去除Fs情况下得到的结果:
图3
由于时域矩形波信号中T=0.4,按理论分析,其频域波形的第⼀个零点数值应该为1/T=2.5,显然,在归⼀化横轴不乘以Fs/2的情况下得到的频谱第⼀零点不是2.5。信号的频域时域不对应。
对于为什么将求得的傅⾥叶变换结果除以N,我感觉百度上⼀位⽹友的回答挺合适,摘抄在下⾯:
因为傅⽴叶变换之后的结果虽然长度和原来数据⼀样
但是前半部分和后半部分结果是共轭对称的
如果只考虑幅度的画,前后两半是关于中⼼对称的
真正有意义的就是0到采样频率⼀半的数据,⽽后半和前半的信息是⼀样的
所以就只取⽤结果的1到N/2,也就是前⼀半的数据
根据变换前后能量相等,原来信号时域上的能量积分和后来信号频域上的积分应该相等
由于只取了⼀半,所以频域的结果能量的积分就会减少⼀半
为了拟补这减少的⼀半,将半信号的幅度根据对称加到前半
所以先取前⼀半的信号,然后在幅度上乘以2,也就拟补了截取⼀半损失的能量
最后,还要将信号除以N的原因是,傅⽴叶变换是个积分变换
写成数学形式的话,是 f(x)dx的积分,实际上函数和⾃变量微分量dx乘积的积分
⽽我们⽤离散信号去计算的时候,只是信号的求和没有乘上x的增量
你可以想像,同样⼀个信号,如果⼀个⽤采样频率Fs采样,得到N的数据
⼀个⽤2*Fs频率采样,就会得到2*N点数据
对着两个信号做离散的傅⽴叶变换,
采样频率⾼,数据点多的信号得到的数值就会⽐采样频率低数据点少的信号⼤⼀倍
为了修正这个问题,所以最终结果除以N
实际上,就是加⼊信号的总时间长度是1,那么N个点,每个点的采样间隔就是1/N
刚才说的计算积分的时候应该乘以积分间隔
所以最后的傅⽴叶变换结果就要乘以1/N,也就是除以N
综上,最后频谱取前⼀半,乘以2,再除以N
他的回答中是对于单边频谱⽽⾔的,因此,得到的傅⾥叶变换结果不仅要除以信号长度L,也即上⽂中的N,同时还要乘以2,MATLAB
help中给出的例⼦中作出的正是单边频谱,所以他将得到的傅⽴叶变换结果除以L,⽽乘以2实在作图时体现的,
plot(f,2*abs(Y(1:NFFT/2+1)))中将频谱信号的幅值乘以了2。
图1中得到的频谱图形是将频谱函数的幅值没有除以L,图3中得到的图像是将频谱函数的结果除以了0.4L,这⾥之所以除以0.4L⽽不直接除以L,是因为尽管定义的信号长度是L,但后⾯的部分全部为0,等效于不存在。⽽幅值不为0的点只有前0.4L,因此这⾥除以0.4L,对⽐理论可知结果是正确的,符合信号频域时域的对应关系。⽽倘若讲频谱信号的幅值除以L,得到的结果如下所⽰,显然是不正确的。
起初,我对信号的频谱幅值除以L
还是挺有疑问的,因为进⾏的是NFFT点的离散傅⽴叶变换,为什么不除以NFFT,反⽽除以原始信号的长度。后来经过思索后发现,此时虽然进⾏离散傅⽴叶变换的点数增加了,但增加的是在原始序列后补充的数值为0的点,类似于前⾯信号频谱幅值应该除以0.4L⽽不是L⼀样,序列后补充的零点是等效于不存在的。况且信号总的能量也并没有增加,因此仍应该除以L。
同时在作出图形时,频域图形中,MATLAB会显⽰0—Fs作为⼀个周期。正常的时域信号的频谱可能是⽆限延长的,但MATLAB仅仅取0—Fs之间的部分。⽽0-Fs
之间的部分具体包含多少个点是有进⾏离散傅⽴叶变换的点数确定的。倘若进⾏的的是NFFT点的离散傅⽴叶变换,那么0—Fs之间就包括NFFT个点。但不能认为包含NFFT个数据点就应该将得到的离散傅⽴叶变换结果除以NFFT,因为离散傅⽴叶变换是对整个时域信号进⾏的,以连续信号为例,时域信号在进⾏离散傅⽴叶变换时进⾏积分变换的是对时域信号进⾏的积分变换,等效到离散信号,时域信
号包含多少个离散点就应该除以多少。⽽NFFT的改变影响的是0-Fs之间显⽰的零点个数。
结论:
得到的矩形波的频谱与理想⽅波函数的理论频谱稍有不同,理性⽅波的理论频谱其在频域⽆限伸展,⽽在MATLAB中只画出了其在0到Fs的部分或者是0到Fs的部分。由于将⽆限延伸的频谱曲线画出来是不现实的,因此,对于双边频谱,MATLAB只画出了在0到Fs的部分;对于单边频谱,MATLAB只画出了在0到Fs/2的部分。同时,在⽤MATLAB对信号进⾏傅⾥叶变换时,应将归⼀化的频域横轴乘以Fs或者
Fs/2,这样得到的频域时域信号才是完美对应的。
附图3的代码:
function Accumulation()
%关于时域⽅波函数的仿真
Fs=100; % Sampling frequency
T=1/Fs; % Sample time
L=100; % Length of signal
t=(0:L-1)*T; % Time vectormatlab求傅里叶变换
%绘制占空⽐为40%的⽅波图形
y1=zeros(1,L); y1(1:0.4*L)=1;
subplot(121); plot(t,y1);
title('占空⽐为40%的矩形波');
xlabel('Time'); ylabel('Amplitude');
%作出对应的频谱图形
NFFT=1024; Y1=fft(y1,NFFT)/(0.4*L);
f1=Fs/2*linspace(0,1,NFFT/2+1);
subplot(122); plot(f1,2*abs(Y1(1:NFFT/2+1))); title('占空⽐为40%的矩形波对应的频谱函数'); xlabel('Frequency (Hz)'); ylabel('|Y(f)|');
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论