hanning_win.m

采样系统实时信号的频谱分析MATLAB程序

% Hanning加窗插值算法
% cgz  2008.11.19

clear;
% 原始数据：直流:0V; 基波：49.5Hz,100V,10deg; HR2:0.5V,40deg;....
hr0=0;f1=50.1;
hr(1)=25*sqrt(2);deg(1)=10;
hr(2)=0;deg(2)=0;
hr(3)=1.755*sqrt(2);deg(3)=40;
hr(4)=0;deg(4)=0;
hr(5)=0.885*sqrt(2);deg(5)=70;
hr(6)=0;deg(6)=0;
hr(7)=1.125;deg(7)=110;
M=7;f=[1:M]*f1; %设定频率

% 采样：N=2048, fs=10kHz.
fs=10000;
N=2048;     % 约10个周期
T=1/fs;
n=[0:N-1];t=n*T;
x=zeros(size(t));
for k=1:M
    x=x+hr(k)*cos(2*pi*f(k)*t+deg(k)*pi/180);
end

% 分析：
w=0.5-0.5*cos(2*pi*n/N);    % hanning window。
Xk=fft(x.*w);
amp=abs(Xk(1:N/2))/N*2;         %幅频
pha=angle(Xk(1:N/2))/pi*180;    %相频
for k=1:N/2
    if(amp(k)<0.01) pha(k)=0;  %当谐波<10mV时，其相位=0
    end
    if(pha(k)<0) pha(k)=pha(k)+360;%调整到0-360度
    end
end
fmin=fs/N;
xaxis=fmin*n(1:N/2);    %横坐标为Hz
kx=round([1:M]*50/fmin);%各次谐波对应的下标(从0开始)
for m=1:M
    km(m)=searchpeaks(amp,kx(m)+1);  %km为谱峰(从1开始)
    if(amp(km(m)+1)<amp(km(m)-1))
        km(m)=km(m)-1;
    end
    beta(m)=amp(km(m)+1)./amp(km(m));
    delta(m)=(2*beta(m)-1)./(1+beta(m));
end
fx=(km-1+delta)*fmin;     %估计频率
hrx=amp(km)*2.*pi.*delta.*(1-delta.*delta)./sin(pi*delta);  %估计幅度
degx=pha(km)-delta.*180/N*(N-1); %估计相位
degx=mod(degx,360); %调整到0-360度
efx=(fx-f)./f*100;          %频率误差
ehr=(hrx-hr)./hr*100;       %幅度误差
edeg=(degx-deg);   %相位误差  

% 结果输出：
subplot(2,2,1);%画出采样序列
plot(t,x);
hold on;
plot(t,x.*w,'r');   %加窗波形
hold off;
xlabel('x(k)');
title('采样序列');
subplot(2,2,2);%画出FFT分析结果
stem(xaxis,amp,'.r');
xlabel('normlized magnitude(*N/2)');
title('FFT分析结果(幅频)');
subplot(2,2,4);
stem(xaxis,pha,'.r');
xlabel('phase(deg)');
title('FFT分析结果(相频)');
subplot(2,2,3);
plot(ehr);
title('幅度误差(%)');

%文本输出
disp '    分析结果的数值见result.txt文件.';
fid=fopen('result.txt','w');
fprintf(fid,'原始数据：f1=%6.1fHz, N=%.f,  fs=%.f \r\n\r\n',f1,N,fs);
fprintf(fid,'谐波次数      1      2      3      4      5      6     7\r\n');
fprintf(fid,'设定频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',f);
fprintf(fid,'估计频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',fx);
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',efx);

fprintf(fid,'设定幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hr);
fprintf(fid,'估计幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hrx);
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',ehr);

fprintf(fid,'设定相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',deg);
fprintf(fid,'估计相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',degx);
fprintf(fid,'误差(度)  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n\r\n',edeg);

%其他数据
fprintf(fid,'谱峰位置理论值：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',[1:M]*f1/fmin);
fprintf(fid,'谱峰位置估计值：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',km-1+delta);
fprintf(fid,'误差（%%）：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',((km-1+delta)-[1:M]*f1/fmin)./([1:M]*f1/fmin)*100);
fprintf(fid,'delta     ：\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',delta);

fclose(fid);