非均匀采样信号的DFT频谱分析实现

guodongliang

对非均匀采样的信号序列 s，和对应的采样时间序列 t，求其DFT，要求频域仍为均匀频点。
用的公式引用汪安民论文中的非均匀采样公式：X（f)=x(t(n))exp(-j*2*pi*f*t(n))*(t(n)-t(n-1));
我自己编写的结果不对，还不如直接进行DFT。
% example
clc
clear
fs=1000;                               % sampling rate  
t=0:1/fs:1;                            % t is equi-spaced
t=t+rand(1,size(t,2))*0.005;% t is not equi-space
f1=50;
f2=100;
f3=0; %200;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);% s is non-uniform sampled sequence
subplot(2,1,1)
plot(t,s)
subplot(2,1,2)
Xk = nudft_gdl(s,t);                        % use nudft-gdl
%Xk = dft(s,length(s));                   % use dft
n=0:length(s)/2;
f = fs*n/length(s);
plot(f,abs(Xk(1:1:length(f))));
title('Non-Uniform DFT')
xlabel('Frequency / Hz')

函数1：离散傅里叶变换DFT的Matlab程序为：
function [Xk]=dft(xn,N)
% compute the DFT coefficients
% [Xk]=dft(xn,N)
% Xk =    DFT coefficients ,0 <= n <=  N-1
% xn = N -length data sequence
% N = length of DFT
%
n=[0:1:N-1];       % raw vector of n
k=[0:1:N-1];       % raw vector of k
WN=exp(-j*2*pi/N); % Wn factor
nk=n'*k;           % matrix
WNnk=WN.^nk;       % DFT matrix
Xk = xn*WNnk;      % raw vector of DFT coefficient

函数2：自己编写的非均匀采样信号的离散傅里叶变换NUDFT的Matlab程序为：
function [Xk]=nudft_gdl(x,t)
% guodl's  compute the dft of non-uniform sampling data
%
% x = non uniform data
% t = time sequency
%
N = length(x);
dt = diff(t);
n1 = [1:N-1];          %t;% % n1=t and k1=n1: no ;%  n1=t and k1=[0:1:N-1]: no  ;%  n1=[0:1:N-1]and k1=t no
k1 = [1:N-1];               %n1;% k的raw vector  % only n1 and k1 are all [0:1:N-1] the result is good
WN=exp(-j*2*pi/(N-1));   % Wn factor
nk1=n1'*k1;           %
WNnk1=WN.^nk1;         % DFT matrix
Xk1 = x(1:N-1)*WNnk1;      %  
Xk = Xk1.*dt;        % NUDFT
我认为n1和k1这两个矢量有问题，应该一个体现均匀的频率格点，一个体现变化时间差。
但仿真发现除非n1 =[1:N-1]；k1=[1:N-1]，其它都不行，如n1=dt;或者n1 = t(1:N-1)

guodongliang

哈哈，我自己搞定了

jeer630

很想知道楼主是怎么搞定的……

massana

回复 1 # guodongliang 的帖子

你好！对非均匀采样我用自己编的程序和你写的程序试过之后，总存在计算出的频率是原信号频率两倍的情况，请问您发现了么？