FFT在功率谱密度计算中的应用
FFT和功率谱估计

FFT和功率谱估计1.用Fourier变换求取信号的功率谱---周期图法clf;Fs=1000;N=256;Nfft=256;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值,并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('周期图 N=256');grid on;Fs=1000;N=1024;Nfft=1024;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值,并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,2),plot(f,Pxx);%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('周期图 N=256');grid on;2.用Fourier变换求取信号的功率谱---分段周期图法%思想:把信号分为重叠或不重叠的小段,对每小段信号序列进行功率谱估计,然后取平均值作为整个序列的功率谱clf;Fs=1000;N=1024;Nsec=256;%数据的长度和FFT所用的数据长度n=0:N-1;t=n/Fs;%采用的时间序列randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第二段功率谱Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第三段功率谱Pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4/4);%Fourier振幅谱平方的平均值,并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('平均周期图(无重叠) N=4*256');grid on;%运用信号重叠分段估计功率谱Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱Pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;%第二段功率谱Pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第三段功率谱Pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;%第四段功率谱Pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第四段功率谱Pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;%第四段功率谱Pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7/7);%Fourier振幅谱平方的平均值,并转化为dBf=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线xlabel('频率/Hz');ylabel('功率谱/dB');title('平均周期图(重叠1/2) N=1024');grid on;3.用Fourier变换求取信号的功率谱---welch方法%思想:welch法采用信号重叠分段,加窗函数和FFT算法等计算一个信号序列的自功率谱(PSD)和两个信号序列的互功率谱(CSD),采用MATLAB自%带的函数psdclf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;window=hanning(256);noverlap=128;dflag='none';randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);Pxx=psd(xn,Nfft,Fs,window,noverlap,dflag);f=(0:Nfft/2)*Fs/Nfft;plot(f,10*log10(Pxx));xlabel('频率/Hz');ylabel('功率谱/dB');title('PSD--Welch方法');grid on;4.功率谱估计----多窗口法(multitaper method ,MTM法)%思想:利用多个正交窗口获得各自独立的近似功率谱估计,综合这些得到一个序列的功率谱估计;相对于普通的周期图有更大的自由度;MTM法采用一个参数:时间带%宽积NW,这个参数用以定义计算功率谱所用窗的数目为2*NW-1,NW越大,时间域分辨率越高而频率分辨率越低,使得功率谱估计的波动减小;随着NW 的增大%,每次估计中谱泄露增多,总功率谱估计的偏差增大clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);[Pxx1,f]=pmtm(xn,4,Nfft,Fs); %此处有问题subplot(2,1,1),plot(f,10*log10(Pxx1));xlabel('频率/Hz');ylabel('功率谱/dB');title('多窗口法(MTM)NW=4');grid on;[Pxx,f]=pmtm(xn,2,Nfft,Fs);subplot(2,1,2),plot(f,10*log10(Pxx));xlabel('频率/Hz');ylabel('功率谱/dB');title('多窗口法(MTM)NW=2');grid on;5.功率谱估计----最大熵法(maxmum entmpy method,MEM法)%思想:假定随机序列为平稳高斯过程利用已知的自相关序列rxx(0),rxx(1),rxx(2)...rxx(p)为基础,外推自相关序列rxx(p+1),rxx(p+2)...保证信息熵最大clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;window=hanning(256);randn('state',0);xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);[Pxx1,f]=pmem(xn,14,Nfft,Fs); %此处有问题subplot(2,1,1),plot(f,10*log10(Pxx1));xlabel('频率/Hz');ylabel('功率谱/dB');title('最大熵法(MEM)Order=14');grid on;%采用Welch方法估计功率谱noverlap=128;dflag='none';subplot(2,1,2)psd(xn,Nfft,Fs,window,noverlap,dflag);xlabel('频率/Hz');ylabel('功率谱/dB');title('Welch方法估计功率谱');grid on;6.功率谱估计----多信号分类法(multiple signal classification,music法)%注:适用于白白噪声中的多正弦波频率估计%思想:将数据自相关矩阵看成是由信号自相关矩阵和噪声自相关矩阵两部分组成,求他们的矩阵特征值向量clf;Fs=1000;N=1024;Nfft=256;n=0:N-1;t=n/Fs;randn('state',0);xn=sin(2*pi*100*t)+2*sin(2*pi*200*t)+randn(1,N);pmusic(xn,[7,1.1],Nfft,Fs,32,16);xlabel('频率/KHz');ylabel('功率谱/dB'); title('Welch方法估计功率谱');grid on;。
功率谱分析及其运用简答题

功率谱分析及其运用简答题一、功率谱分析的基本原理功率谱分析的基本思想是将一个连续时间的信号转换为频域上的离散信号,然后对这些离散信号进行傅里叶变换,得到其频谱表示。
频谱表示中的每个峰值代表了一个特定的频率分量,而每个峰值的高度则代表了该频率分量的强度。
通过对频谱表示进行加权平均,可以得到原始信号的能量分布情况。
二、功率谱分析的应用场景1.通信系统:在无线通信系统中,功率谱分析可以用来检测干扰信号或者识别出合法的通信信号。
通过比较接收到的信号与已知的噪声信号之间的功率谱差异,可以判断出是否存在干扰。
此外,功率谱分析还可以用来估计信道容量和误码率等重要参数。
2.音频处理:在音频处理中,功率谱分析可以用来提取音乐中的基音和谐波等信息。
通过对音乐信号进行快速傅里叶变换(FFT),可以得到其频谱表示,然后再通过滤波器等算法提取出所需的信息。
3.雷达系统:在雷达系统中,功率谱分析可以用来检测目标反射回来的信号。
通过对反射回来的信号进行功率谱分析,可以确定目标的位置、速度和形状等信息。
三、实际运用举例下面以一个简单的示例来说明功率谱分析的实际运用过程。
假设我们有一个包含多个正弦波成分的信号x(t),我们需要将其分解成若干个简单的正弦波成分y(i),并计算每个成分的振幅和频率。
具体步骤如下:1.对信号x(t)进行快速傅里叶变换(FFT),得到其频域表示f (k)。
2.对频域表示f(k)进行平滑处理,以减少高频噪声的影响。
常用的平滑方法包括均值滤波和中值滤波等。
3.对平滑后的频域表示f(k)进行平方运算,得到其功率谱密度ρ(f)。
4.根据需要,可以选择不同的窗函数对ρ(f)进行加窗处理,以减少频谱泄漏等问题。
常见的窗函数包括汉宁窗、汉明窗和矩形窗等。
5.最后,根据ρf)的大小和位置等信息,可以确定原始信号中包含的各个正弦波成分以及它们的振幅和频率等特征。
cpm信号功率谱密度的快速计算方法

标题:CPM信号功率谱密度的快速计算方法摘要:本文主要探讨了连续相位调制(CPM)信号功率谱密度的快速计算方法。
首先介绍了CPM信号的基本原理和特点,然后详细讨论了传统计算方法的局限性和缺点。
提出了一种基于快速傅里叶变换(FFT)的高效计算方法,并对该方法进行了数学推导和仿真验证。
对比了传统方法和提出的快速计算方法的优劣,并给出了进一步研究的展望。
关键词:CPM信号,功率谱密度,快速计算,傅里叶变换正文:一、引言连续相位调制(CPM)是一种在通信系统中广泛应用的调制方式,其可以有效地提高系统的频谱利用率和抗多径干扰能力。
在CPM系统中,信号的功率谱密度对系统性能有着重要的影响。
对CPM信号功率谱密度的快速计算方法具有重要的理论意义和实际应用价值。
二、CPM信号的基本原理和特点CPM信号是一种通过改变载波相位来传输信息的调制方式。
与传统的调频调相调制方式相比,CPM信号具有带宽效率高、抗多径干扰能力强等优点。
其调制框图如下图所示:(这里可以插入CPM信号的调制框图)三、传统计算方法的局限性和缺点传统的计算CPM信号功率谱密度的方法一般是通过信号的自相关函数来推导。
这种方法在理论上是可行的,但在实际计算中存在着以下几个方面的局限性和缺点:1. 计算复杂度高:传统方法需要进行复杂的积分计算和级数求和,计算复杂度较高;2. 计算速度慢:由于计算复杂度高,传统方法的计算速度较慢,不适用于实时性要求较高的场景;3. 数值稳定性差:由于计算中涉及积分和级数求和,容易出现数值稳定性差的情况,影响计算结果的准确性。
传统计算方法存在着计算复杂度高、计算速度慢和数值稳定性差等局限性和缺点。
四、基于FFT的高效计算方法为了克服传统计算方法的局限性和缺点,本文提出了一种基于快速傅里叶变换(FFT)的高效计算方法。
该方法的基本思想是将CPM信号的功率谱密度计算问题转化为频域上的复杂信号处理问题,利用FFT 算法对信号进行高效处理,从而实现功率谱密度的快速计算。
matlab计算功率谱密度时fft点数、窗函数的作用

matlab计算功率谱密度时fft点数、窗函数的作用
在MATLAB中,计算功率谱密度通常涉及到傅里叶变换。
FFT(快速傅里叶变换)是实现这一目的的关键工具。
以下是FFT点数和窗函数在计算功率谱密度中的作用:
1. **FFT点数**:
* FFT点数决定了频谱的分辨率。
更多的点数意味着更精细的频率分辨率,但也会增加计算复杂性。
* 在进行傅里叶变换时,你需要选择一个FFT点数。
例如,如果你选择了256个点,那么你可以分辨到256Hz的频率,这是因为一个FFT周期为256个样本点。
* 在大多数应用中,为了获得准确的频率分辨率,FFT长度应该是信号长度的2倍或更高的倍数。
2. **窗函数**:
* 窗函数在傅里叶变换之前应用于信号,有助于减少频谱泄漏。
频谱泄漏是由于信号的突然开始和结束在频域产生的高频成分导致的。
* 窗函数可以在信号上加一个窗口,该窗口在开始和结束处逐渐变为零,从而减少信号边缘的影响。
* 不同的窗函数(如汉明窗、汉宁窗、矩形窗等)有不同的特性,包括边缘斜率、主瓣宽度和旁瓣高度。
选择合适的窗函数可以减少旁瓣,从而提高频率分辨率。
使用窗函数和正确的FFT点数可以更准确地估计信号的功率谱密度,
尤其是在处理非平稳信号时。
在MATLAB中,你可以使用`fft`函数进行傅里叶变换,并使用pwelch函数(功率谱估计函数)来估计功率谱密度,其中你可以指定窗函数和FFT点数。
数字信号处理技术在信号检测中的应用

数字信号处理技术在信号检测中的应用概述数字信号处理技术在现代通信系统中扮演着至关重要的角色。
在数字信号处理中,我们将信号抽样、量化/编码、数字滤波、FFT等方法应用到信号处理中。
数字信号处理技术可以广泛应用于语音、图像处理、雷达、通信等领域,提高了信号处理的精度和速度。
在本文中,我们将探讨数字信号处理技术在信号检测中的应用。
信号检测技术及其分类信号检测是信号处理领域的一个重要分支,用于确定通过信道传输的信号是否存在。
信号检测技术可以分为基于时间域和基于频域的方法。
时间域的方法针对时间序列信号进行操作,例如信号的差分、平均等,在信号的功率谱密度不明显时适用。
频域方法则将信号转换为频域上的函数,例如将信号通过FFT算法转换为时频图,在信号的功率谱密度较明显时适用。
数字滤波数字滤波是数字信号处理中最常见的技术之一。
数字滤波可以分为时域滤波和频域滤波。
时域滤波针对时间序列信号,在时域上进行卷积运算,例如低通滤波器和高通滤波器,可以用于去除噪声、平滑信号和保留信号的一定频率分量;频域滤波器则将信号转换到频域上,通过乘上某些频率分量来去除部分信号分量,例如带通滤波器、陷波和带阻滤波器,可以用于去除干扰噪声和选择特定频率分量。
FFT快速傅里叶变换(FFT)是将时域信号转换成频域信号的一种方法。
FFT算法不仅可以用来分析频域上的信号,还可以用来压缩数据和进行频域上的滤波处理。
在信号检测中,FFT技术可以用来分析信号在频域上的特征,例如特定频率分量的能量。
通过对信号进行FFT变换,可以更准确地分析信号特征和区分噪声信号和有效信号。
数字信号处理在信号检测中的应用数字信号处理技术可以广泛应用于信号检测领域。
下面我们将探讨数字信号处理技术在雷达、语音处理和信号处理中的应用。
雷达信号处理雷达信号处理是将雷达回波信号转换为原始数据、成像和目标识别的关键技术。
利用数字信号处理技术,我们可以提高雷达回波信号的分辨率和精度。
例如,可以使用带宽滤波器来过滤杂乱的回波信号,并通过FFT算法提取目标的频域特征。
fft计算的功率谱

fft计算的功率谱
快速傅里叶变换(FFT)是一种用于计算离散傅里叶变换(DFT)和其逆变换的算法。
在信号处理中,FFT常被用于将信号从时域转换到频域,从而方便我们分析信号的频率成分。
功率谱密度(Power Spectral Density, PSD)描述的是信号或者时间序列的频率内容。
对于一个信号,其功率谱给出了信号在各个频率上的功率分布。
计算功率谱的一个常用方法就是使用FFT。
以下是使用FFT计算功率谱的基本步骤:
1. 采集信号:首先,你需要有一个时域信号。
这个信号可以是一段时间内的声音、电压等物理量的测量值。
2. 应用窗函数:为了减小频谱泄漏,通常在信号上应用一个窗函数,如汉宁窗或海明窗。
3. 执行FFT:对加窗后的信号执行快速傅里叶变换,得到频域表示。
4. 计算功率谱:对FFT的结果取模平方,然后除以信号长度N(或者乘以2除以N,这取决于你的FFT实现和是否需要归一化),得到功率谱。
如果信号是单
边带(只考虑正频率),那么还需要乘以2来保留总功率。
数学上,这个过程可以表示为:
X[k]=FFT(x[n])X[k] = FFT(x[n])X[k]=FFT(x[n])
Pxx[k]=∣X[k]∣2NPxx[k] = \frac{|X[k]|^2}{N}Pxx[k]=N∣X[k]∣2
其中,xxx[n]x[n]x[n] 是时域信号,XXX[k]X[k]X[k] 是其频域表示,Pxx[k]Pxx[k]Pxx[k] 是在频率kkk 处的功率谱密度。
注意:这只是计算功率谱的一种基本方法,实际应用中可能还需要考虑其他因素,如窗函数的选择、重叠处理、平滑等。
matlab 功率谱密度计算

matlab 功率谱密度计算
Matlab是一种功能强大的计算机软件,能够完成许多数学和工程计算任务。
其中一个常见的应用是计算功率谱密度。
功率谱密度是一个信号在频域上的能量分布,通常用于频域分析和信号处理。
在 Matlab 中,可以使用 fft 函数将一个信号从时域转换到频域。
转换后,可以使用 abs 函数计算信号频谱的幅度,然后将其平方以计算信号的功率谱密度。
最后,可以使用 plot 函数将功率谱密度显示为频率的函数。
以下是一个简单的 Matlab 代码示例,用于计算并显示信号的功率谱密度:
% 导入信号数据
data = importdata('signal.dat');
% 计算信号的频域表示
freq_domain = fft(data);
% 计算信号的功率谱密度
psd = abs(freq_domain).^2;
% 显示功率谱密度图
plot(psd);
xlabel('Frequency');
ylabel('Power Spectral Density');
这是一个基本的示例,可以根据需要进行修改和扩展。
Matlab 提供了许多功能,可用于处理和分析各种类型的信号数据。
matlab中计算功率谱的4种方法

在MATLAB中,计算功率谱是信号处理和频谱分析中的重要任务。
功率谱可以帮助我们了解信号中不同频率成分的能量分布情况,对于理解信号特性和进行频谱分析都是至关重要的。
在MATLAB中,有多种方法可以用来计算功率谱,在本文中,我将介绍并比较其中的四种常用方法。
第一种方法是使用MATLAB中的`periodogram`函数。
`periodogram`函数可以直接计算信号的功率谱密度(PSD),它采用傅里叶变换的方法,将信号从时域转换到频域,并计算功率谱密度。
这种方法简单直接,适用于对功率谱快速估计的情况。
在使用`periodogram`函数时,我们可以指定窗函数和重叠比例等参数,来对功率谱的估计进行优化。
第二种方法是使用`pwelch`函数。
`pwelch`函数也可以用来计算信号的功率谱密度,它采用Welch方法,通过对信号进行分段,然后对每个段进行傅里叶变换,并对结果进行平均来估计功率谱密度。
Welch 方法可以减小估计的方差,得到更平滑和可靠的功率谱估计结果。
在使用`pwelch`函数时,同样可以指定窗函数和重叠比例等参数来优化估计结果。
第三种方法是使用`fft`函数和自行计算功率谱。
通过对信号进行傅里叶变换得到频谱,然后对频谱的幅度进行平方运算,即可得到功率谱。
这种方法的好处是灵活性高,可以根据具体需求对傅里叶变换和求平方的结果进行后续处理,比如进行平滑或滤波操作。
但是需要注意的是,自行计算功率谱需要对信号处理和频谱分析有较深的理解。
第四种方法是使用`cpsd`函数。
`cpsd`函数可以用来计算信号之间的交叉功率谱密度,适用于多信号系统中不同信号之间的频谱分析。
交叉功率谱密度可以帮助我们理解不同信号之间频率成分的相关性和影响程度,对于系统建模和故障诊断都是非常有帮助的。
MATLAB提供了多种方法来计算功率谱,每种方法都有其适用的场景和优势。
在具体应用中,我们可以根据信号特性和分析需求来选择合适的方法。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
F F T在功率谱密度计算中的应用集团标准化工作小组 #Q8QGGQT-GX8G08Q8-GNQGJ8-MHHGN#FFT在功率谱密度计算中的应用一、FFT算法理论依据和编程思想FFT算法的基本思想:考察DFT与IDFT的运算发现,利用以下两个特性可减少运算量:Ⅰ)系数是一个周期函数,它的周期性和对称性可利用来改进运算,提高计算效率。
如:因此利用这些周期性和对称性,DFT运算中有些项可合并;Ⅱ)利用W N nk的周期性和对称性,把长度为N点的大点数的DFT运算分解为若干个小点数的DFT。
因为DFT的计算量正比于N2,N小计算量也就小。
FFT算法正是基于这样的基本思想发展起来的。
它有多种形式,下面是按时间抽取的FFT(N点DFT运算的分解)先从一个特殊情况开始,假定N是2的整数次方,N=2M,M:正整数1.将N点的DFT分解为两个N/2点的DFT:首先将序列x(n)分解为两组,一组为偶数项,一组为奇数项r=0,1,…,N/2-1将DFT运算也相应分为两组:其中X1(k)和X2(k)分别是x1(r)和x2(r)的N/2点DFT。
可见,一个N点的DFT可以分解为两个N/2点的DFT,这两个N/2点的DFT再按照上面(1)式合成为一个N点DFT,注意到,X1(k),X2(k)有N/2个点,即k=0,1,…,N/2-1,由(1)式得到X(k)只有N/2点,而实际上X(k)有N个点,即k=0,1,…,N-1,要用X1(k),X2(k)表示全部X(K)值,还必须应用系数w的周期性和对称性。
(k)的(N/2)~N-1点表示:由X(k)= X1(k)+w kNX2(k), k=0,1,2,…,N/2-1得:(2a),因为,且同样。
k对称性:。
考虑到WN故(2b)(2a)式表示了X(k)前半部分k=0~N/2-1时的组成方式,(2b)式则表示了后半部分k=N/2~N-1时的组成方式。
这两式所表示的运算过程可用一个称作蝶形的信号流图来表示。
3.蝶形信号流图:如图1(a)所示,图中左面两支为输入,中间以一个小圆圈表示加、减运算,右上支为相加输出,右下支为相减输出,如果在某一支路上信号需要进行乘法运算,则在该支路上标以箭头,并将相乘的系数标在简头边,这样(2a),(2b)所表示的运算,可用图1(b)所表示的“蝶形结”来表示。
采用这种表示法,可将以上以讨论的分解过程用计算流图来表示。
图所示为N=23=8的例子。
通过这样分解以后,每一个N/2点DFT只需要图 N点DFT分解为2个N/2点DFT(N=8)(N/2)2= N2/4次复数乘法运算,两个N/2点的DFT需要2(N/2)2= N2/2 次复乘,再加上将两个N/2点DFT合成为N点DFT时,在蝶形结前的N/2次复乘,共需要(N/2)2+N/2≈ N2/2次复乘,由此可见,经过这样的分解处理,运算量差不多节省了一倍。
4.将N/2点的DFT分解为两个N/4点的DFT:既然这样分解是有效的,由于N=2M,N/2仍然是偶数,因此可对两个N/2点的DFT再分别作进一步分解,例如对x1(r)和x2(r)可以再按其偶数部分及奇数部分分解为两个N/4点的DFT,既然这样分解是有效的,由于N=2M,N/2仍然是偶数,因此可对两个N/2点的DFT再分别作进一步分解,例如对x1(r)和x2(r)可以再按其偶数部分及奇数部分分解为两个N/4点的DFT,l=0,1,…,N/4-1而同样X2(k)也可这样分解,并且将系数统一为,这样一个8点DFT就可分解为四个2点的DFT,如图所示。
图 N点DFT分解为4个N/4点的DFT(N=8)个点DFT的表示:最后剩下的是2点DFT,它可以用一个蝶形结表示,例如,x(0),x(4)所组成的2点DFT就可表示式:这样,一个8点的完整的按时间抽取运算的流图如图所示。
由于这样的方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的,所以称为“按时间抽取法”或“时间抽取法”。
6.时间抽取法FFT运算特点:(1)蝶形运算对任何一2的整数幂N=2M,总是可以通过M次分解最后完全成为2点的DFT运算。
这样的M次分解,就构成从x(n)到X(k)的M级运算过程。
从上面的流图可看到,每一级运算都由N/2个蝶形运算构成。
因此每一级运算都需要(N/2次复乘和N次复加(每个结作加、减各一次),这样,经过时间抽取后M级运算总共需要的运算:复乘复加N·M=Nlog2N实际运算量与这个数字稍有出入,因为W这几个系数实际上都不用乘法运算,因此在上面N=8的例子中,实际上只有两个系数WN 1及WN3是需要乘法运算的。
用时间抽取法所需的计算量,不论是复乘还是复加都与Nlog2N成正比,而直接运算时则与N2成正比。
例N=2048,N2=4194304,(N/2)log2N=11264,N2/[(N/2)log2N]=倍。
FFT显然要比直接法快得多。
(2)原位计算当数据输入到存储器中以后,每一级运算的结果仍然储存在同一组存储器中,直到最后输出,中间无需其它存储器,这叫原位计算。
例如,N=8的FFT运算,输入x(0),x(4),x(2),x(6)…,x(7)可分别存入A(1),A(2),…,A(8)这9个存储单元中,在第一级运算中,首先是存储单元A(1),A(2)中x(0),x(4)进入蝶形运算,x(0),x(4)输入运算器后,其数值不再需要保存,因此蝶形运算的结果可仍然送回存储单元A(1),A(2)中保存,然后A(3),A(4)中x(2),x(6)再进入蝶形运算,其结果再送回A(3),A(4),一直到算完A(7),A(8),则完成了第一级运算过程。
第二级运算仍可采用这种原位的方式,但是进入蝶形结的组合关系不同,首先进入蝶形结的是A(1)、A(3)存储单元中的数据,运算结果仍可送回A(1)、A(3)保存,然后进入蝶形结的是A(2)、A(4)…,依此类推,每一级运算均可在原位进行,这种原位运算结构可节省存储单元,降低设备成本,还可节省找地址的时间。
(3)序数重排对按时间抽取FFT的原位运算结构,当运算完毕时,这种结构存储单元A(1)、A(2),…,A(8)中正好顺序存放着X(0),X(1),X(2),…,X(7),因此可直接按顺序输出,但这种原位运算的输入x(n)却不能按这种自然顺序存入存储单元中,而是按X(0),X(4),X(2),X(6),…,X(7)的顺序存入存储单元,这种顺序看起来相当杂乱,然而它也是有规律的。
当用二进制表示这个顺序时,它正好是“码位倒置”的顺序。
例如,原来的自然顺序应是x(1)的地方,现在放着x(4),用二进制码表示这一规律时,则是在x(0 0 1)处放着x(1 0 0),x(0 1 1)处放着x(1 1 0)。
即将自然顺序的二进制码位倒置过来,第一位码变成最末位码,这样倒置以后的顺序正是输入所需要的顺序。
下表列出N=8时按码位倒置规律所得的顺序,其结果与按时间抽取FFT流图中的输入顺序是一致的。
表一码位倒置顺序自然顺序二进码表示码位倒置码位倒置顺序00000000在实际运算中,一般直接将输入数据x (n )按码位倒置的顺序排好输入很不方便,总是先按自然顺序输入存储单元,然后再通过变址运算将自然顺序的存储转换成码位倒置顺序的存储,然后进行FFT 的原位计算。
目前有许多通用DSP 芯片支持这种码位倒置的寻址功能。
(4)蝶形类型随迭代次数成倍增加观察8点FFT 的三次迭代运算第一级迭代,只有一种类型的蝶形运算系数W 08第二级迭代,有二种类型的蝶形运算系数W 08、W 28,参加运算的两个数据点间隔为2。
第三级迭代,有四类蝶形运算系数W 08、W 18、W 28、W 38,参加运算的两个数据点间隔为4。
所以,每次迭代的蝶形类型比上一次蝶代增加一倍,数据点间隔也增大一倍。
7.功率谱密度的计算根据相关定理与维纳-辛钦关系式可得随机信号序列x(n)的功率谱密度 21()lim ()x N s k X k N→∞= (1) 其估计值21()()x S k X k N∧=(2) 如果观察到序列x(n)的N 个值,即x(0), x(1),…x(N-1),就可以通过FFT 直接求得X(k),再按式(2)求得S x (k),其计算过程如图所示。
1.程序框图(1) 运算主程序框图 图周期图法计算功率谱密度流程N=1<<M;jiaodu=(-2*PI)/N;////进行fft计算: ///fft();int zengzhi;double t1r,t2r,t1i,t2i;double jd; ///计算的角度int wz=1<<M;for( i=0;i<N;i++){xr[i]=x[i];xi[i]=0;}for(int i1=M-1;i1>=0;i1--){wz=wz/2;kuaishu=1<<(i1);geshu=N/kuaishu;zengzhi=kuaishu;for(int i2=0;i2<kuaishu;i2++){for(int i3=0;i3<(geshu/2);i3++){jd=jiaodu*zengzhi*i3;t1r=xr[geshu*i2+i3];t1i=xi[geshu*i2+i3];t2r=xr[geshu*i2+i3+geshu/2];t2i=xi[geshu*i2+i3+geshu/2];xr[geshu*i2+i3]=t1r+t2r*cos(jd)-t2i*sin(jd);xi[geshu*i2+i3]=t1i+t2i*cos(jd)+t2r*sin(jd);xr[geshu*i2+i3+geshu/2]=t1r-t2r*cos(jd)+t2i*sin(jd);xi[geshu*i2+i3+geshu/2]=t1i-t2i*cos(jd)-t2r*sin(jd);}}}for(i=0;i<512;i++){shibu[i]=xr[i];xubu[i]=xi[i];}}四、计算实例采样函数:x(t)=A*cos(2*PI*f1*t)+B*sin(2*PI*f2*t)采样点数:N=512计算结果:(打开程序运行时保存的文件操作过程:1.打开生成的exe文件(),则会出现以下窗口:2.点击软件界面的“参数导入”按钮,就可以计算完成了,生成文件保存在所在文件夹。
计算结果为:X[0]=+0 j X[1]=+ j X[2]=+ jX[3]=+ j X[4]=+ j X[5]=+ jX[6]=+ j X[7]=+ j X[8]=+ jX[9]=+ j X[10]=+ j X[11]=+ jX[12]=+ j X[13]=+ j X[14]=+ jX[15]=+ j X[16]=+ j X[17]=+ jX[18]=+ j X[19]=+ j X[20]=+ jX[21]=+ j X[22]=+ j X[23]=+ jX[24]=+ j X[25]=+ j X[26]=+ jX[27]=+ j X[28]=+ j X[29]=+ jX[30]=+ j X[31]=+ j X[32]=+ jX[33]=+ j X[34]=+ X[35]=+ jX[36]=+ j X[37]=+ X[38]=+ jX[39]=+ j X[40]=+ j X[41]=+ jX[42]=+ j X[43]=+ X[44]=+ jX[45]=+ j X[46]=+ j X[47]=+ jX[48]=+ j X[49]=+ X[50]=+ jX[51]=+ j X[52]=+ X[53]=+ jX[54]=+ j X[55]=+ X[56]=+ jX[57]=+ j X[58]=+ X[59]=+ jX[60]=+ j X[61]=+ X[62]=+ jX[63]=+ j X[64]=+ j X[65]=+ jX[66]=+ j X[67]=+ j X[68]=+ jX[69]=+ j X[70]=+ j X[71]=+ jX[72]=+ j X[73]=+ X[74]=+ jX[75]=+ j X[76]=+ X[77]=+ jX[78]=+ j X[79]=+ j X[80]=+ jX[81]=+ j X[82]=+ X[83]=+ jX[84]=+ j X[85]=+ X[86]=+ jX[90]=+ j X[91]=+ X[92]=+ j X[93]=+ j X[94]=+ X[95]=+ j X[96]=+ j X[97]=+ X[98]=+ j X[99]=+ j X[100]=+ X[101]=+ j X[102]=+ j X[103]=+ X[104]=+ j X[105]=+ j X[106]=+ X[107]=+ j X[108]=+ j X[109]=+ X[110]=+ j X[111]=+ j X[112]=+ X[113]=+ j X[114]=+ j X[115]=+ X[116]=+ j X[117]=+ j X[118]=+ X[119]=+ j X[120]=+ j X[121]=+ X[122]=+ j X[123]=+ j X[124]=+ X[125]=+ j X[126]=+ j X[127]=+ X[128]=+ j X[129]=+ j X[130]=+ X[131]=+ j X[132]=+ j X[133]=+ X[134]=+ j X[135]=+ j X[136]=+ X[137]=+ j X[138]=+ j X[139]=+ X[140]=+ j X[141]=+ j X[142]=+ X[143]=+ j X[144]=+ j X[145]=+ X[146]=+ j X[147]=+ j X[148]=+ X[149]=+ j X[150]=+ j X[151]=+ X[152]=+ j X[153]=+ j X[154]=+ X[155]=+ j X[156]=+ j X[157]=+ X[158]=+ j X[159]=+ j X[160]=+ X[161]=+ j X[162]=+ j X[163]=+ X[164]=+ j X[165]=+ j X[166]=+ X[167]=+ j X[168]=+ j X[169]=+ X[170]=+ j X[171]=+ j X[172]=+ X[173]=+ j X[174]=+ j X[175]=+ X[176]=+ j X[177]=+ j X[178]=+ X[179]=+ j X[180]=+ j X[181]=+ X[182]=+ j X[183]=+ j X[184]=+ X[185]=+ j X[186]=+ j X[187]=+ X[188]=+ j X[189]=+ j X[190]=+ X[191]=+ jX[195]=+ j X[196]=+ X[197]=+ j X[198]=+ j X[199]=+ X[200]=+ j X[201]=+ j X[202]=+ X[203]=+ j X[204]=+ j X[205]=+ X[206]=+ j X[207]=+ j X[208]=+ X[209]=+ j X[210]=+ j X[211]=+ X[212]=+ j X[213]=+ j X[214]=+ X[215]=+ j X[216]=+ j X[217]=+ X[218]=+ j X[219]=+ j X[220]=+ X[221]=+ j X[222]=+ j X[223]=+ X[224]=+ j X[225]=+ j X[226]=+ X[227]=+ j X[228]=+ j X[229]=+ X[230]=+ j X[231]=+ j X[232]=+ X[233]=+ j X[234]=+ j X[235]=+ X[236]=+ j X[237]=+ j X[238]=+ X[239]=+ j X[240]=+ j X[241]=+ X[242]=+ j X[243]=+ j X[244]=+ X[245]=+ j X[246]=+ j X[247]=+ X[248]=+ j X[249]=+ j X[250]=+ X[251]=+ j X[252]=+ j X[253]=+ X[254]=+ j X[255]=+ j X[256]=+0j X[257]=+ j X[258]=+ j X[259]=+ X[260]=+ j X[261]=+ j X[262]=+ X[263]=+ j X[264]=+ j X[265]=+ X[266]=+ j X[267]=+ j X[268]=+ X[269]=+ j X[270]=+ j X[271]=+ X[272]=+ j X[273]=+ j X[274]=+ X[275]=+ j X[276]=+ j X[277]=+ X[278]=+ j X[279]=+ j X[280]=+ X[281]=+ j X[282]=+ j X[283]=+ X[284]=+ j X[285]=+ j X[286]=+ X[287]=+ j X[288]=+ j X[289]=+ X[290]=+ j X[291]=+ j X[292]=+ X[293]=+ jX[294]=+ j X[295]=+ X[296]=+ jX[297]=+ j X[298]=+ X[299]=+X[300]=+ j X[301]=+ j X[302]=+ j X[303]=+ j X[304]=+ j X[305]=+ j X[306]=+ j X[307]=+ j X[308]=+ j X[309]=+ j X[310]=+ j X[311]=+ j X[312]=+ j X[313]=+ j X[314]=+ j X[315]=+ j X[316]=+ j X[317]=+ j X[318]=+ j X[319]=+ j X[320]=+ j X[321]=+ j X[322]=+ j X[323]=+ j X[324]=+ j X[325]=+ j X[326]=+ j X[327]=+ j X[328]=+ j X[329]=+ j X[330]=+ j X[331]=+ j X[332]=+ j X[333]=+ j X[334]=+ j X[335]=+ j X[336]=+ j X[337]=+ j X[338]=+ j X[339]=+ j X[340]=+ j X[341]=+ j X[342]=+ j X[343]=+ j X[344]=+ j X[345]=+ j X[346]=+ j X[347]=+ j X[348]=+ j X[349]=+ j X[350]=+ j X[351]=+ j X[352]=+ j X[353]=+ j X[354]=+ j X[355]=+ j X[356]=+ j X[357]=+ j X[358]=+ j X[359]=+ j X[360]=+ j X[361]=+ j X[362]=+ j X[363]=+ j X[364]=+ j X[365]=+ j X[366]=+ j X[367]=+ j X[368]=+ j X[369]=+ j X[370]=+ j X[371]=+ j X[372]=+ j X[373]=+ j X[374]=+ j X[375]=+ j X[376]=+ j X[377]=+ j X[378]=+ j X[379]=+ j X[380]=+ j X[381]=+ j X[382]=+ j X[383]=+ j X[384]=+ j X[385]=+ j X[386]=+ j X[387]=+ j X[388]=+ j X[389]=+ j X[390]=+ j X[391]=+ j X[392]=+ j X[393]=+ j X[394]=+ j X[395]=+ j X[396]=+ j X[397]=+ j X[398]=+ j X[399]=+ j X[400]=+ j X[401]=+ jX[402]=+ j X[403]=+ j X[404]=+ j X[405]=+ j X[406]=+ j X[407]=+ j X[408]=+ j X[409]=+ j X[410]=+ j X[411]=+ j X[412]=+ j X[413]=+ j X[414]=+ j X[415]=+ j X[416]=+ j X[417]=+ j X[418]=+ j X[419]=+ j X[420]=+ j X[421]=+ j X[422]=+ j X[423]=+ j X[424]=+ j X[425]=+ j X[426]=+ j X[427]=+ j X[428]=+ j X[429]=+ j X[430]=+ j X[431]=+ j X[432]=+ j X[433]=+ j X[434]=+ j X[435]=+ j X[436]=+ j X[437]=+ j X[438]=+ j X[439]=+ j X[440]=+ j X[441]=+ j X[442]=+ j X[443]=+ j X[444]=+ j X[445]=+ j X[446]=+ j X[447]=+ j X[448]=+ j X[449]=+ j X[450]=+ j X[451]=+ j X[452]=+ j X[453]=+ j X[454]=+ j X[455]=+ j X[456]=+ j X[457]=+ j X[458]=+ j X[459]=+ j X[460]=+ j X[461]=+ j X[462]=+ j X[463]=+ j X[464]=+ j X[465]=+ j X[466]=+ j X[467]=+ j X[468]=+ j X[469]=+ j X[470]=+ j X[471]=+ j X[472]=+ j X[473]=+ j X[474]=+ j X[475]=+ j X[476]=+ j X[477]=+ j X[478]=+ j X[479]=+ j X[480]=+ j X[481]=+ j X[482]=+ j X[483]=+ j X[484]=+ j X[485]=+ j X[486]=+ j X[487]=+ j X[488]=+ j X[489]=+ j X[490]=+ j X[491]=+ j X[492]=+ j X[493]=+ j X[494]=+ j X[495]=+ j X[496]=+ j X[497]=+ j X[498]=+ j X[499]=+ j X[500]=+ j X[501]=+ j X[502]=+ j X[503]=+ j X[504]=+ j X[505]=+ j X[506]=+ jX[507]=+ j X[508]=+ j X[509]=+ jX[510]=+ j X[511]=+ j X[512]=3+ j3.分析计算结果点击软件界面的“功率谱密度图”按钮,显示如下图所示:。