zfft.m

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

fs=1000;
N=4096;
D=50;
M=200;

t=(0:N*D+2*M)/fs;
x=4*sin(2*pi*166.35*t+0.8*sin(2*pi*0.65*t)+pi/6)+0.6*randn(1,N*D+2*M+1);

xf=fft(x,N);
xf=abs(xf(1:N/2))/N*2;

subplot(121);
plot((0:N/2-1)*fs/N,xf);
axis([164 170 0 4])
grid on
xlabel('时间(t)');ylabel('幅值(g)')
set(gca,'XTick',[164 166 168]);

fe=167;
L=1;  %不做平均

k=1:M;                          
w=0.5+0.5*cos(pi*k/M);          %Hanning窗

fl=max(fe-fs/(4*D),-fs/2.2);
fh=min(fe+fs/(4*D),fs/2.2);

yf=D*fl;                       %移频量
df=fs/D/N;
f=fl:df:fl+(N/2-1)*df;
xz=zeros(1,N/2);
wl=2*pi*fl/fs;
wh=2*pi*fh/fs;
hr(1)=(wl-wh)/pi;
hr(2:M+1)=(sin(wl*k)-sin(wh*k))./(pi*k).*w;
hi(1)=0;
hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;

k=0:N-1;
w=0.5-0.5*cos(2*pi*k/N);

for i=1:L
    for k=1:N
        kk=(k-1)*D+M+(i-1)*N;
        xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
        xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
    end
    xzt=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf/fs);
    xzt=xzt.*w;
    xzt=xzt-sum(xzt)/N;
    xzt=fft(xzt);
    xz=xz+(abs(xzt(1:N/2))/N*2).^2;
end
xz=(xz/L).^0.5;
subplot(122);
plot((0:N/2-1)*fs/N,xf,f,xz);
axis([163.5 169 0 4])
grid on
xlabel('频率(Hz)');ylabel('幅值(g)')