加窗相位差频率matlab,求助：关于相位差校正，加卷积窗函数

本帖最后由 芒果樱桃 于 2012-4-8 09:05 编辑

希望实现功能是通过加自卷积窗函数，获得正确的相位、频率和幅值

根据丁康的相位差校正，试着编了一个程序，将信号分成两段，然后分别加窗，通过相位差校正法求出信号信息，

不知程序这样写行不行，新手求解答

%相位差校正

clc

clear all

fs=1024;

N=1024;

M=N/2;

L=512;

t1=0:(N-1);  %序列一段长度

t2=0:(N+L-1);  %序列二段长度(本来是从L到L+N)

df=fs/N;%频谱分辨率

%信号函数

ph1=0.1;%相位1

ph3=0;%相位3

f1=50.2;

f3=99.8;

x=0.2+6*sin(2*pi*f1*t1/N+ph1*pi/180)+sin(2*pi*f3*t1/N+ph3*pi/180);

x2=0.2+6*sin(2*pi*f1*t2/N+0.1*pi/180)+sin(2*pi*f3*t2/N+ph3*pi/180);

%plot(linspace(0,pi,N),angle(fft(x,N)));

%原始信号函数的输出

X=fft(x);  %做FFT变换

Ayy=abs(X)/(N/2);  %取模，换算成实际的幅度

Ayy(1)=Ayy(1)/2;  %直流分量是(1),上句话多乘了2，所以这句要除以2

f=(1:N-1)*fs/N;  %换算成实际的频率值

for k=1:N/2

ph(k)=phase(X(k))*180/pi+90; %计算相位，换算为角度

end;

%加blackmanharris自卷积窗

window=blackmanharris(M)';

win=conv(window,window);  %blackmanharris窗函数的时域相乘，即为频域卷积

win1=ifft(fft(win,N));

win2=[zeros(1,length(x2)-N),win];%对窗函数补零，序列2段本来应该也是长度为N

win2=ifft(fft(win2,N+L));

X1=fft(x.*win1);  %信号1段加窗

A1=abs(X1);

Ph1=angle(X1)*180/pi+90;  %求相位

Ph1=[Ph1,zeros(1,L)];

X2=fft(x2.*win2);

A2=abs(X2)/(N/2);

Ph2=angle(X2)*180/pi+90;

dph=Ph2-Ph1;

delta=dph-2*pi*L/N;

delta=mod(delta,2*pi);

for k=1:round(f/df);%小数四舍五入，这里返回值是整数,k为谱线号

if  delta(k)>pi %delta为两序列谱峰线间的相位差，(0~2*pi)

delta(k)=delta(k)-2*pi;

elseif delta(k)

delta(k)=delta(k)+2*pi;

end

end

dfy=(dph-2*pi*L*k/N)/(2*pi*L/N);

fconr=(k+dfy)*df;

sita=atan(imag(X)/real(X))+dfy;