Recursive Time-Frequency Reassignment

Recursive Time-Frequency Reassignment published in IEEE Transactions on Signal Processing

Today I proudly present the original manuscript from my 2009 IEEE Transactions on Signal Processing publication. It can be downloaded by clicking the article thumbnail. I also share a corresponding MATLAB implementation in the code snippets below.

The functions rspec() and rrspec() in rspec.m and rrspec.m computes the recursive Short-Time Fourier Transform (STFT) and the reassigned version, respectively. A usage example is included in the file usage_ex.m: it shows how to apply the functions on a linear chirp in the 40 to 80 Hz bandwidth.

 

Usage Example (usage_ex.m):

%
% Usage example for the functions rspec() and rrspec()
% This script can be found online at http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2008
%

fs = 1000;
t=0:0.001:1;
s=chirp(t,40,.5,60);
s = [zeros(1,100) s zeros(1,100)];

rtfr = rspec(s, fs, 4, 1:.5:100);
rrtfr = rrspec(s, fs, 4, 1:.5:100);

subplot(2,1,1), imagesc(rtfr), subplot(2,1,2), imagesc(rrtfr);

 

Recursive STFT spectrogram (rspec.m):

function rtfr = rspec(x,fs,k,f);
% rtfr = rtfr(x,fs,k,f) computes the recursive STFT spectrogram with 
% the constants defined
% x - input signal
% fs - sampling frequency
% k - order
% f - frequency vector
%
% Equations referred to are found in the journal paper "Recursive
% Time-Frequency Reassignment" available in IEEE's Transactions on Signal
% Processing. This script can be found online at http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2009
% Updated December, 2012.
%

if k < 2 || k > 5
    sprintf('k must be in the range 2 <= k <= 5')
    rtfr = 0;
    return;
end;

T = 1/fs;
nf = length(f);
N = length(x);
rtfr = zeros(nf, N);

Omega = f(2)-f(1); % frequency spacing

% Balance time-frequency resolution in number of samples, eq. 16
sigma_p = (sqrt(Omega)*factorial(k-1))/(sqrt(2*pi*T)*(k-1)^(k-1)*exp(-(k-1)));
coef = cell(10,1);

for j = 1:nf;
    omega_p = 2*pi*f(j);
    p = -sigma_p + i*omega_p; % eq. 6
    a = exp(p*T);
 
    % Coefficients from eq. 14, k = 2...5
    % Numerator coefficients in coef{i}, i = 2n+1, n = 1,2,3,4
    % Demoninator coefficients in coef{i}, i = 2n, n = 1,2,3,4,5

    coef{3} = sigma_p^2*T^2*[0 a];
    coef{4} = [1 -2*a a.^2];
    coef{5} = sigma_p^3*T^3*[0 1/2*a 1/2*a.^2];
    coef{6} = [1 -3*a 3*a.^2 -a.^3];
    coef{7} = sigma_p^4*T^4*[0 1/6*a 2/3*a.^2 1/6*a.^3];
    coef{8} = [1 -4*a 6*a.^2 -4*a.^3 a.^4];
    coef{9} = sigma_p^5*T^5*[0 1/24*a 11/24*a.^2 11/24*a.^3 1/24*a.^4];
    coef{10}= [1 -5*a 10*a.^2 -10*a.^3 5*a.^4 -a.^5];
  
    rtfr(nf-j+1,:) = filter(coef{2*k-1},coef{2*k},x); % Eq. 16
end;

rtfr = abs(rtfr);

 

Recursive Reassigned STFT spectrogram (rrspec.m):

function [rrtfr coef dcoef tcoef] = rrspec(x,fs,k,f);
% rrtfr = rrspec(x,fs,k,f) computes the recursive reassigned STFT spectrogram with 
% the constants defined
% x - input signal
% fs - sampling frequency
% k - order
% f - frequency vector
%
% Equations referred to are found in the journal paper "Recursive 
% Time-Frequency Reassignment" available in IEEE's Transactions on Signal 
% Processing. This script can be found online at http://8-void.com
%
% Author: Geir Kjetil Nilsen (geir.kjetil.nilsen@gmail.com) 2009
% Updated December, 2012.
%
if k < 2 || k > 4
 sprintf('k must be in the range 2 <= k <= 4')
 rrtfr = 0;
 return;
end;

T = 1/fs;
nf = length(f);
N = length(x);
w = zeros(nf, N);
w_d = zeros(nf, N);
w_t = zeros(nf, N);
rrtfr = zeros(nf,N);

Omega = f(2)-f(1);

% Balance time-frequency resolution in number of samples, eq. 16
sigma_p = (sqrt(Omega)*factorial(k-1))/(sqrt(2*pi*T)*(k-1)^(k-1)*exp(-(k-1))); 

coef = cell(2,1);
dcoef = cell(2,1);
tcoef = cell(2,1);

for j = 1:nf;
    omega_p = 2*pi*f(j);
    p = -sigma_p + i*omega_p; % eq. 6
    a = exp(p*T);

    % Coefficients from eq. 14, k = 2...5
    
    % Numerator coefficients in coef{i}, i = 2n+1, n = 1,2,3,4
    % Demoninator coefficients in coef{i}, i = 2n, n = 2,3,4,5

    coef{3} = sigma_p^2*T^2*[0 a];
    coef{4} = [1 -2*a a.^2];
    coef{5} = sigma_p^3*T^3*[0 1/2*a 1/2*a.^2];
    coef{6} = [1 -3*a 3*a.^2 -a.^3];
    coef{7} = sigma_p^4*T^4*[0 1/6*a 2/3*a.^2 1/6*a.^3];
    coef{8} = [1 -4*a 6*a.^2 -4*a.^3 a.^4];
    coef{9} = sigma_p^5*T^5*[0 1/24*a 11/24*a.^2 11/24*a.^3 1/24*a.^4];
    coef{10}= [1 -5*a 10*a.^2 -10*a.^3 5*a.^4 -a.^5];

    % Coefficients from eq. 23, k = 2...4
    dcoef{3} = sigma_p.^2*T*[1 -a*(1 - p*T)];
    dcoef{4} = coef{4}; 
    dcoef{5} = sigma_p.^3*T.^2*[0 a/2*(2 + p*T) -a.^2/2*(2 - p*T)];
    dcoef{6} = coef{6}; 
    dcoef{7} = sigma_p.^4*T.^3*[0 a/6*(3 + p*T) a.^2*2/3*p*T -a.^3/6*(3 - p*T)];
    dcoef{8} = coef{8};

    % Coefficients from eq. 24, k = 2...4
    tcoef{3} = 2/sigma_p*coef{5}; 
    tcoef{4} = coef{6}; 
    tcoef{5} = 3/sigma_p*coef{7}; 
    tcoef{6} = coef{8}; 
    tcoef{7} = 4/sigma_p*coef{9}; 
    tcoef{8} = coef{10};

    w(j,:) = filter( coef{2*k-1}, coef{2*k},x); % Eq. 16
    w_d(j,:) = filter(dcoef{2*k-1}, dcoef{2*k},x); % Freq. reass. coeffs with eq. 16
    w_t(j,:) = filter(tcoef{2*k-1}, tcoef{2*k},x); % Time reass. coeffs with eq. 16
end;

rtfr = abs(w);
step = f(2)-f(1); 
f0 = f(1);

% eq. 27
for n = 1:N;
    for f = 1:nf;
    that = n - round(1/T*real(w_t(f,n)./w(f,n))); % eq. 26
    that = min(max(that, 1),N); % edge
 
    fhat = round(1/step*imag(w_d(f,n)./w(f,n))/2.0/pi - 1/step*(f0)); % eq. 25
    fhat = mod(fhat, nf) + 1; % edge
 
    rrtfr(nf - fhat + 1,that) = rrtfr(nf - fhat + 1,that) + rtfr(f,n); 
    end;
 end;