全相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

全相FFT

相关推荐