zoomfft.m

这里是两个matlab开发的频率细化分析算法

%细化FFT
clear
close all hidden
format long

fni=input('细化FFT处理-输入数据文件名 :','s');
fid=fopen(fni,'r');         %以只读方式打开数据文件
sf=fscanf(fid,'%f',1);      %读入采样频率值
fi=fscanf(fid,'%f',1);      %最小细化截至频率
np=fscanf(fid,'%d',1);      %放大倍数
nfft=fscanf(fid,'%d',1);    %FFT长度
x=fscanf(fid,'%f',[1,inf]); %读入时程数据存成行向量
status=fclose(fid);         %关闭数据文件

nt=length(x);
%最大细化截至频率(
fa=fi+0.5*sf/np;
%取大于n且最接近n的2的整数次方为FFT的长度
nf=2^nextpow2(nt);
%确定细化带宽的数据长度
na=round(0.5*nf/np+1);

%频移
%建立一个按1递增的向量
n=0:nt-1;
%确定单位旋转因子向量(b=2*pi*ft*t,ft为细化频带中心频率，即ft=(fi+fa)/2)
b=n*pi*(fi+fa)/sf;
%乘以单位旋转因子进行频移
y=x.*exp(-i*b);

%频移后的低通滤波                     %        滤波和下采样
%FFT变换                             %      c=decimate(y,np); 
b=fft(y,nf);                         %
%正频率带通内的元素赋值                %  左边的语句可用上面的一句代替！
a(1:na)=b(1:na);                     %
%负频率带通内的元素赋值                %
a(nf-na+1:nf)=b(nf-na+1:nf);         %
%IFFT变换                            %
b=ifft(a,nf);                        %
                                     %
%重采样                              %
c=b(1:np:nt);                        %

%进行细化FFT变换
a=fft(c,nfft)*2/nfft;
%变换结果重新排序
                %FFT对称是有条件的，即实数序列的FFT将是实部偶对称，虚部奇对称。
                %但在exzfft_m中，为了细化，把输入序列x乘以旋转因子进行频移：
                %b=n*pi* (fi+fa)/fs;     %乘单位旋转因子进行频移
                %y=x.*exp(-i*b);
                %y便是一个复数序列了，即使经过滤波和下采样，依然是一个复数序列，y的FFT便是复数序列的FFT，FFT后不再满足对称的性质了。
                %同样在细化以后，采样频率为sf/np，信号复调 制后0频率为原信号的 (fi+fa)/2，信号一样满足采样定理，
                %信号的最高频率应小于sf/np/2，对应于原信号便是在fi~fa之间，即对应于y2的第四部分和第一部分。
y2=zeros(1,nfft/2);
%排列负频率段的数据
y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);
%排列正频率段的数据
y2(nfft/4+1:nfft/2)=a(1:nfft/4);
n=0:(nfft/2-1);
%定义细化FFT频率向量
f2=fi+n*2*(fa-fi)/nfft;

%对输入数据做FFT用来变换比较
%FFT变换
y1=fft(x,nfft)*2/nfft;

f1=n*sf/nfft;
%定义与细化FFT相同的频率范围
ni=round(fi*nfft/sf+1);
na=round(fa*nfft/sf+1);

%绘制输入曲线图形
subplot(211);
t=0:1/sf:(nt-1)/sf;
nn=1:1024;
plot(t(nn),x(nn))
xlabel('时间(s)');ylabel('幅值');grid on;
%在相同的频率范围下绘制曲线图
subplot(212);
nn=ni:na;
plot(f1(nn),abs(y1(nn)/length(y1)),':',f2,abs(y2)/length(y2));
xlabel('频率(Hz)');ylabel('幅值');
legend('未细化','细化');grid on;