%% 比值法频谱校正——第五章
%% 产生仿真信号
clear
N=256;
fs=256;
t=0:1/fs:(N-1)/fs;
y=cos(2*pi*50.4*t+pi/4);
a=0.5; %加n窗
T0=N/(fs);
t=0:1/fs:(N-1)/fs;
w=a-(1-a)*cos(2*pi*t/T0);
Y1=y.*w;
Kt=2;
%% 做信号的fft
%Y=fft(Y);
W=exp(-j*2*pi/N); %自编fft算法实现
M=ceil(log(N)/log(2));
x=0:N-1; %位码倒序
b=dec2bin(x,M);
c=fliplr(b);
x=bin2dec(c);
Y=zeros(1,N); %将信号值按位码倒序重新排列
for n=0:N-1
Y(n+1)=Y1(x(n+1)+1);
end
for i=1:M %快速傅里叶运算
for j=0:2^(i-1)-1
for n=0:2^i:N-1
Y(n+j+1)=Y(n+j+1)+Y(n+2^(i-1)+j+1)*W^(N*j/2^i); %实现蝶形运算
Y(n+2^(i-1)+j+1)=Y(n+j+1)-2*Y(n+2^(i-1)+j+1)*W^(N*j/2^i);
end
end
end
n=0:N-1;
y(n+1)=2*abs(Y(n+1))/N; %得幅值
y2=y;
for n=0:N-1 %得相位角
if y(n+1)<0.1
Y(n+1)=0;
else Y(n+1)=Y(n+1);
end
end
Xr=real(Y);
Xi=imag(Y);
for n=0:N-1
if Xr(n+1)>0
y1(n+1)=atan(Xi(n+1)/Xr(n+1))/pi*180;
else
if Xi(n+1)<0
y1(n+1)=atan(Xi(n+1)/Xr(n+1))/pi*180-180;
else
y1(n+1)=atan(Xi(n+1)/Xr(n+1))/pi*180+180;
end
end
end
f=0:fs/N:fs-fs/N;
figure(3)
subplot(2,2,1);stem(f,y2) %输出幅值谱
title('比值频谱校正法')
xlabel('频率(HZ)-校正前')
ylabel('幅值')
axis([0 100 0 1.2]);
subplot(2,2,3);stem(f,y1) %输出相位谱
xlabel('频率(HZ)-校正前')
ylabel('相位（°）')
axis([0 100 -180 180]);
%% 采用比值法校正
a=0.1; %寻找幅值的峰值,a为阈值
for n=0:N-1
if y(n+1)<+a
y(n+1)=a;
else y(n+1)=y(n+1);
end
end
[yk k]=findpeaks(y);
yk=yk(1);
k=k(1)-1;
if y(k)>=y(k+2) %找次峰值
yk1=y(k);
v=yk1/yk;
df1=(v-2)/(v+1);
f01=k-(df1+1);
else
yk1=y(k+2);
v=yk1/yk;
df1=(v-2)/(v+1);
f01=k+(df1+1);
end
f0=f01*fs/N; %校正频率
df1=k-f01;
df=df1*fs/N;
A=2*pi*df1*yk*(1-df1^2)/sin(pi*df1); %校正幅值
S=y1(k+1)+180*df1; %校正相位
f=0-df:fs/N:fs-fs/N-df;
y2(k+1)=A;
y1(k+1)=S;
figure(3)
subplot(2,2,2);stem(f,y2) %输出校正幅值谱
xlabel('频率(HZ)-校正后')
ylabel('幅值')
axis([0 100 0 1.2]);
subplot(2,2,4);stem(f,y1) %输出校正相位谱
xlabel('频率(HZ)-校正后')
ylabel('相位（°）')
axis([0 100 -180 180]);
