注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

Second generation Vold-Kalman Order Filtering  

2012-06-16 14:03:57|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

The Vold-Kalman Filter, introduced by H?vard Vold and Jan Leuridan in 1993, is able to extract non-stationary periodic components from a signal using a known frequency vector [1]. Being formulated in a least-squares sense, it can be solved as a sparse linear system. Similar to the Kalman filter, VKF minimises a cost function based on a structural equation and a data equation.

This submission implements the second generation VKF with the ability to extract multiple orders at the same time, with energy spreading in case of crossing orders

 

 

function varargout = vkf(y,fs,f,p,bw)
%VKF 2nd Generation Vold-Kalman Order Filtering.
%   x = VKF(y,fs,f) extracts the order with frequency vector f from signal
%   y with samplerate fs, using a 2-pole filter with a -3dB bandwidth of
%   1 percent of the sample rate. The output is a single waveform x.
%
%   [...] = VKF(y,fs,f,p) uses a p-order filter (typically between 1 or 4).
%   Every order increases the roll-off by -40dB per decade. By specifying
%   additional lower-order coefficients, zero boundary conditions are
%   added. For instance: p = [2 0 1] applies 2nd order filtering and
%   forces the envelope amplitude and its first derivative to zero at t_1
%   and t_N.
%
%   [...] = VKF(y,fs,f,p,bw) uses a bandwidth in Hertz specified by bw. If
%   bw is a scalar, a constant bandwidth is used; if bw is a vector with
%   the same length as y, a time-varying instantaneous bandwidth is
%   realised.
%
%   X = VKF(y,fs,F,...) with [N,K] = size(F), performs simultaneous
%   extraction of K orders with frequency vectors [f_1,...,f_K] in array
%   F. In case of crossing orders, this method tries to reveal the correct
%   order amplitudes. The output is an array of K waveforms [x_1,...,x_K].
%
%   [a,c] = VKF(...) returns the complex envelope(s) a and phasor(s) c,
%   such that the order waveform(s) can be reconstructed by x = real(a.*c).
%
%   [a,c,r] = VKF(...) ouputs an additional selectivity vector r used
%   to realise the bandwidth given by bw.
%
%   Note: Filter orders > 4 usually result in ill conditioning and should
%   be avoided. The filter bandwidth determination was implemented for
%   arbitrary order but was not verified for orders higher than 3.
%
%   Demo:
%   Calling VKF without arguments shows a small demonstration of multi-
%   order filtering with two crossing orders in the presence of white
%   noise. Note that the demo uses the spectrogram function from the Signal
%   Processing Toolbox.
%
%   Example:
%       fs = 4000;
%       T = 5;
%       dt = 1/fs;
%       t = (0:dt:(T-dt))';
%       N = numel(t);
%
%       % Instationary component
%       A1 = [linspace(0.5,1,floor(N/2)) linspace(1,0.5,ceil(N/2))]';
%       f1 = 0.2*fs - 0.1*fs*cos(pi*t/T);
%       phi1 = 2*pi*cumsum(f1)*dt;
%       y1 = A1.*cos(phi1);
%
%       % Stationary component
%       A2 = ones(N,1);
%       f2 = 0.2*fs*ones(N,1);
%       phi2 = 2*pi*cumsum(f2)*dt;
%       y2 = A2.*sin(phi2);
%
%       % White noise
%       e = 2*rand(size(y1));
%
%       % Mixed signal
%       y = y1 + y2 + e;
%
%       % Perform VKF on periodic components
%       p = 2;
%       bw = 1;
%       [a,c] = vkf(y,fs,[f1 f2],p,bw);
%       x = real(a.*c);
%
%       % Reveal white noise
%       w = y-sum(x,2);
%
%       % Plot
%       figure('color','white')
%       subplot(2,1,1)
%       spectrogram(y,round(fs/16),[],[],fs)
%       subplot(2,1,2), hold on
%       plot(t,[A1 A2],'--')
%       plot(t,abs(a))
%       xlabel('Time (s)')
%       ylabel('Amplitude')
%
%       % Playback
%       soundsc(y,fs) % Original signal
%       soundsc(x,fs) % Periodic components (2 channels)
%       soundsc(w,fs) % Remaining white noise
%
%   Written by: Maarten van der Seijs, 2010.
%   Version 1.3, 20 April 2012.
%
%   References:
%   [1] Vold, H. and Leuridan, J. (1993), High resolution order tracking
%       at extreme slew rates, using Kalman tracking filters. Technical
%       Report 931288, Society of Automotive Engineers.
%
%   [2] Tuma, J. (2005), Setting the passband width in the Vold-Kalman
%       order tracking filter. Proceedings of the International Congress
%       on Sound and Vibration (ICSV12), Lisbon, Portugal.


%% Input processing

if (nargin == 0) && (nargout == 0), vkfdemo(), return, end
if nargin < 3, error(generatemsgid('ArgChk'),'Wrong number of input arguments. Use: vkf(y,fs,f).'); end
if nargin < 4, p = 2; end           % Default filter order
if nargin < 5, bw = fs/100; end     % Default bandwidth: 0.01*fs

if size(y,2) > size(y,1), y = transpose(y); end
if size(f,2) > size(f,1), f = transpose(f); end
if size(bw,2) > size(bw,1), bw = transpose(bw); end

% Signal
n_t = length(y);

% Frequency vector(s)
[n_f, n_ord] = size(f);
if n_f == 1
    f = ones(n_t,1)*f;
elseif n_f ~= n_t
    error(generatemsgid('ArgChk'),'The vectors in f should have 1 or %d elements.',n_t)
end

% Turn frequency vectors into phasor vectors
c = exp(2i*pi*cumsum(f,1)/fs);

% Bandwidth vector
n_bw = numel(bw);
if n_bw == 1
    bw = ones(n_t,1)*bw;
elseif n_bw ~= n_t
    error(generatemsgid('ArgChk'),'Vector bw should have 1 or %d elements.',n_t)
end

% Relative bandwidth in radians
phi = pi/fs.*bw;

% Filter order
p_tmp = p;
p = max(p);
p_lo = setdiff(p_tmp,p);


%% Construct filter matrix and bandwidth vector

% Coefficients of p-order difference equation
P_lu = pascal(p+1,1);
coeffs = P_lu(end,:);

% Linear system of difference equations
A = spdiags(ones(n_t-p,1)*coeffs,0:p,n_t-p,n_t);

% Introduce lower order equations to set boundary conditions to zero
pp = numel(p_lo);
A_pre = sparse(pp,n_t);
A_pre(:,1:p) = P_lu(p_lo+1,1:p);
A = [A_pre; A; A_pre(end:-1:1,end:-1:1)];

% Determination of filter r-value
s = 0:p;
sgn = (-1).^s;

% Allocate p+1 linear system
Q = zeros(p+1,p+1);
b = zeros(p+1,1);

% Sum of coefficients
Q(1,:) = ones(p+1,1);
b(1) = 2^(2*p-1);

% System of p diff. equations
for i = 1:p
    Q(i+1,:) = s.^(2*(i-1)) .* sgn;
end

% Solve for q
q = round(Q\b).' .* sgn;

% Calculate r
num = sqrt(2)-1;
den = zeros(size(bw));
for qi = 1:numel(q)
    den = den + 2*q(qi)*cos((qi-1)*phi);
end
r = sqrt(num./den);

% Check potential conditioning of B-matrix
if any(den <= 0) || any(r > sqrt(1/(2*q(1)*eps)))
    error(generatemsgid('IllConditioned'),'Ill-conditioned B-matrix; selectivity bandwidth is too small.')
end

% Generate single-order B matrix
R = spdiags(r,0,n_t,n_t);
B = (A*R)'*(A*R) + speye(n_t);

% Free memory
clear A A_pre R bw phi num den f

%% Construct multi-order matrices

% Construct sparse diagonal part
nn = n_t*n_ord;
diags = -p:p;
diags_B = spdiags(B,diags);
BB_D = spdiags(repmat(diags_B,n_ord,1),diags,nn,nn);

% Prepare sparse upper-diagonal part
bl_U = (n_ord^2 - n_ord)/2;
ii_U = zeros(n_t,bl_U);
jj_U = zeros(n_t,bl_U);
cc_U = complex(zeros(n_t,bl_U),0);
m = 0;
% Upper-diagonal part
for ki = 1:n_ord
    for kj = (ki+1):n_ord
        m = m + 1;
        ii_U(:,m) = (ki-1)*n_t + (1:n_t);
        jj_U(:,m) = (kj-1)*n_t + (1:n_t);
        cc_U(:,m) = conj(c(:,ki)).*c(:,kj);       
    end   
end

% Construct sparse upper-diagonal part
BB_U = sparse(ii_U,jj_U,cc_U,nn,nn);

% Assemble sparse matrix
BB = BB_D + BB_U + BB_U';

% Construct right-hand side
cy = conj(reshape(c,nn,1)).*repmat(y,n_ord,1);

% Free memory
clear diags_B B BB_D BB_U ii_U jj_U cc_U

%% Solve

xx = BB\cy;
x = 2 * reshape(xx,[n_t,n_ord]);

switch nargout
    case 1 % Output real-valued order waveforms
        varargout{1} = real(x.*c);
    case 2 % Output complex-valued envelopes and phasors
        varargout{1} = x;
        varargout{2} = c;
    case 3 % Output additional r-vectors
        varargout{1} = x;
        varargout{2} = c;
        varargout{3} = r;
end
end

function vkfdemo()
fs = 4000;
T = 5;
dt = 1/fs;
t = (0:dt:(T-dt))';
N = numel(t);

% Instationary component
A1 = [linspace(0.5,1,floor(N/2)) linspace(1,0.5,ceil(N/2))]';
f1 = 0.2*fs - 0.1*fs*cos(pi*t/T);
phi1 = 2*pi*cumsum(f1)*dt;
y1 = A1.*cos(phi1);

% Stationary component
A2 = ones(N,1);
f2 = 0.2*fs*ones(N,1);
phi2 = 2*pi*cumsum(f2)*dt;
y2 = A2.*sin(phi2);

% White noise
e = 2*rand(size(y1));

% Mixed signal
y = y1 + y2 + e;

% Perform VKF on periodic components
p = 2;
bw = 1;
[a,c] = vkf(y,fs,[f1 f2],p,bw);
x = real(a.*c);

% Reveal white noise
w = y-sum(x,2);

% Plot
figure('color','white')
subplot(2,1,1)
spectrogram(y,round(fs/16),[],[],fs)
subplot(2,1,2), hold on
plot(t,[A1 A2],'--')
plot(t,abs(a))
xlabel('Time (s)')
ylabel('Amplitude')

% Playback
soundsc(y,fs) % Original signal
soundsc(x,fs) % Periodic components (2 channels)
soundsc(w,fs) % Remaining white noise
end

  评论这张
 
阅读(1018)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016