**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;