全相FFT

作者：桂。

时间：2017-12-02  23:29:48

链接：http://www.cnblogs.com/xingshansi/p/7956491.html 

---

一、相位提取

以正弦信号为例，x = sin(2pi*f*t+pi)，希望提取phi：

思路1：通过Hilbert变化解决

思路2：借助FFT，利用插值思想，估计Phi;

思路3：借助全相FFT（apFFT, all phase FFT）实现。

思路三可提取信号相位，这一点FFT做不到，而相位信息通常可判断相位调制类型，可用于情报的脉内检测。

全相FFT思路：

• 选定N点窗，如hanning
• 窗函数自相关，并归一化
• 对2N-1序列x(n)加窗，
• 将2N-1个点，每间隔N点相加
• FFT实现

二、仿真验证

clc;clear all;close all;

fs = 1e9;
fo =  200e6;
t = 0:1/fs:1023/fs;
tao = 0.3/3e8*sind(45);
SNR = 20;
ch1 = awgn(sin(2*pi*t*fo),SNR) ;
ch2 = awgn(sin(2*pi*t*fo + 2*pi*tao*fo),SNR);
% ch1 = sin(2*pi*t*fo) ;
% ch2 = sin(2*pi*t*fo + 0.5);
pha = angle(hilbert(ch2))-angle(hilbert(ch1));
figure()
subplot 211
plot(t*fs,pha);
subplot 212
plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on;
%FFT提取相位
pha1 = angle(fft(ch1(1:512)).*fft(ch2(1:512)));
figure()
subplot 211
plot(t(1:512)*fs,pha1);
subplot 212
plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on;
%apFFT提取相位
win = hanning(512)';
win1 = conv(win,win);
win1 = win1/sum(win1);
y1 = ch1(1:1023).*win1;
y2 = ch2(1:1023).*win1;
out1 = [0,y1(1:511)]+y1(512:1023);
out2 = [0,y2(1:511)]+y2(512:1023);
pha2 = angle(fft(out1).*conj(fft(out2)));
figure()
subplot 211
plot(t(1:512)*fs,pha2);
subplot 212
plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on;
[~,pos] = max(abs(fft(ch1(1:512))));
[pha2(pos) mean(pha) ;-pi+ 2*pi*tao*fo  2*pi*tao*fo]
theta_est = asind((pha2(pos))/2/pi/fo/0.3*3e8)+90;
abs(theta_est-45)

另外，频谱细化，可借助zoom-FFT：

fs = 2048;
T = 100;
t = 0:1/fs:T;
x = 30 * cos(2*pi*110.*t) + 30 * cos(2*pi*111.45.*t) + 25*cos(2*pi*112.3*t) + 48*cos(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t);
[f, y] = zfft(x, 109, 115, fs);
plot(f, y);


function [f, y] = zfft(x, fi, fa, fs)
% x为采集的数据
% fi为分析的起始频率
% fa为分析的截止频率
% fs为采集数据的采样频率
% f为输出的频率序列
% y为输出的幅值序列(实数）

f0 = (fi + fa) / 2;              %中心频率
N = length(x);                 %数据长度

r = 0:N-1;
b = 2*pi*f0.*r ./ fs;               
x1 = x .* exp(-1j .* b);          %移频

bw = fa - fi;                                       

B = fir1(32, bw / fs);             %滤波 截止频率为0.5bw
x2 = filter(B, 1, x1);               

c = x2(1:floor(fs/bw):N);           %重新采样
N1 = length(c);
f = linspace(fi, fa, N1);
y = abs(fft(c)) ./ N1 * 2;                         
y = circshift(y, [0, floor(N1/2)]);            %将负半轴的幅值移过来
end

　　

参考：Digital Receiver-based Electronic Intelligence System Configuration for the Detection and Identification of Intrapulse Modulated Radar Signals