%?Hanning加窗插值算法
%?cgz??2008.11.19

clear;
%?原始數據：直流:0V;?基波：49.8Hz100V10deg;?HR2:0.5V40deg;....
hr0=0;f1=49.8;
hr（1）=1;deg（1）=pi/12;
hr（2）=2;deg（2）=pi/3;
hr（3）=3;deg（3）=pi/7;
hr（4）=9;deg（4）=pi/6;
hr（5）=4;deg（5）=pi/4;
hr（6）=6;deg（6）=2*pi/3;
hr（7）=8;deg（7）=pi;
M=7;f=[1:M]*f1;?%設定頻率

%?采樣：N=1024?fs=5000Hz.
fs=5000;
N=1024;????
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

%加入高頻噪聲
y1=x+5*cos（2*pi*1000*t）;
%濾波
Y1=fft（y1）;Y1=fftshift（Y1）/1024;
N1=length（y1）;
w1=2*[0:N1-1]./N1-1;f_1=w1/（2*T）;
figure（1）
plot（f_1abs（Y1））;grid;

%濾波器
alphamax=1;alphamin=10;
wp=0.2;ws=0.22;
[N2wh]=buttord（wpwsalphamaxalphamin）;
[ba]=butter（N2wh）;
y2=filter（bay1）;
Y2=fft（y2）;Y2=fftshift（Y2）/512;
figure（2）
plot（f_1abs（Y2））;grid;

%?分析：
x=y2;
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（ampkx（m）+1）;??%km為譜峰（從1開始）
????if（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（degx360）;?%調整到0-360度
efx=（fx-f）./f*100;??????????%頻率誤差
ehr=（hrx-hr）./hr*100;???????%幅度誤差
edeg=（degx-deg）;???%相位誤差??

%?結果輸出：
figure（3）
subplot（221）;%畫出采樣序列
plot（tx）;
hold?on;
plot（tx.*w‘r‘）;???%加窗波形
hold?off;
xlabel（‘x（k）‘）;
title（‘采樣序列‘）;
subplot（222）;%畫出FFT分析結果
stem（xaxisamp‘.r‘）;
xlabel（‘normlized?magnitude（*N/2）‘）;
title（‘FFT分析結果（幅頻）‘）;
subplot（224）;
stem（xaxispha‘.r‘）;
xlabel（‘phase（deg）‘）;
title（‘FFT分析結果（相頻）‘）;
subplot（223）;
plot（ehr）;
title（‘幅度誤差（%）‘）;

%文本輸出
disp?‘????分析結果的數值見result.txt文件.‘;
fid=fopen（‘result.txt‘‘w‘）;
fprintf（fid‘原始數據：f1=%6.1fHz?N=%.f??fs=%.f?\r\n\r\n‘f1Nfs）;
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

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件?????????933??2017-12-20?19:45??result.txt
?????文件?????????422??2008-11-19?19:35??searchpeaks.m
?????文件????????3585??2017-12-20?19:45??hanning_win.m