function demo_bpframemul
%DEMO_BPFRAMEMUL Frame multiplier acting as a time-varying bandpass filter
%
% This demo demonstrates creation and effect of a Gabor multiplier. The
% multiplier performs a time-varying bandpass filtering. The band-pass
% filter with center frequency changing over time is explicitly created
% but it is treated as a "black box" system. The symbol is identified by
% "probing" the system with a white noise and dividing DGT of the output
% by DGT of the input element-wise. The symbol is smoothed out by a 5x5
% median filter.
%
% Figure 1: The symbol of the multiplier.
%
% This figure shows a symbol used in the Gabor multiplier.
%
% Figure 2: Spectroram obtained by re-analysis of the test signal after applying the multiplier
%
% This figure shows a spectrogram of the test signal after applying
% the estimated Gabor multiplier.
%
%
% Url: http://ltfat.github.io/doc/demos/demo_bpframemul.html
% Copyright (C) 2005-2023 Peter L. Soendergaard <peter@sonderport.dk> and others.
% This file is part of LTFAT version 2.6.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% Sampling rate
fs = 44100;
% Input signal
f = gspi;
% Input length
L = numel(gspi);
% Initial center frequency of the band-pass filter
fc = 5000;
% Modulation parameter for the center frequency of the band-pass filter
fvar = 4800;
% Gabor system parameters
win = {'tight','hann',882};
a = 200;
M = 1000;
% The bandpass filter chnges it's center frequency every 2*a samples
Lch = ceil(L/(2*a));
L = Lch*2*a;
% Center frequencies
l = 0:Lch-1;
l = l(:);
fcm = fc+fvar*sin(2*pi*l/Lch*10);
% Parameters of the bandpass filter
G = 40;
Q = 4;
Ha = zeros(Lch,3);
Hb = zeros(Lch,3);
% Design a IIR peaking filters ...
for l=1:Lch
[Ha(l,:),Hb(l,:)]=parpeak(fcm(l),Q,G,fs);
end
% ...and make them band-pass filters.
Ha = Ha/10^(G/20);
% Create a white noise and divide it into input into blocks
ff = randn(L,1);
fblocks = reshape(ff,numel(ff)/Lch,Lch);
% Do a blockwise filtering of the white noise
y = blockfilter(fblocks,Ha,Hb);
% Obtain DGT of the output ...
cy = dgtreal(y(:),win,a,M);
% And the input
cx = dgtreal(fblocks(:),win,a,M);
% Create a symbol by a plain division of DGT of the output and the input.
% (This would be a deconvolution if Fourier multiplier was used)
symbol = cy./cx;
% Smooth out the extremes using a 2D 5x5 median filter
symbol = cmedfilt2(symbol,[5,5]);
% Plot the symbol
figure(1);
plotdgtreal(symbol,a,M,44100,'dynrange',50);
% Apply the symbol to a test signal
f = postpad(f,L);
c = dgtreal(f,win,a,M);
c = c.*symbol;
fhat = idgtreal(c,win,a,M,L);
% Plot spectrogram of the resulting signal
figure(2);
c = dgtreal(fhat,win,a,M);
plotdgtreal(c,a,M,44100,'dynrange',50);
% Filter the input signal by the time-varying band-pass filter directly
y = blockfilter(reshape(f,numel(f)/Lch,Lch),Ha,Hb);
% Export the signals (since we are in a function)
assignin('base','forig', f(:));
assignin('base','fmul', fhat(:));
assignin('base','fdir', y(:));
disp('The original signal can be played by typing: sound(forig,44100);');
disp(['The signal obtained by a direct filtering cen be played by ',...
'(this is what is approximated by the multiplier): sound(fdir,44100);']);
disp(['The signal obtained by applying the Gabor multiplier: sound(fmul,44100);']);
function y = blockfilter(fb,Ha,Hb)
%BLOCKFILTER Block filtering
y = zeros(size(fb));
Lch = size(fb,2);
Z = zeros(2,1);
for l=1:Lch
[y(:,l),Z] = filter(Ha(l,:),Hb(l,:),fb(:,l),Z);
end
function fo = cmedfilt2(f,d)
% CMEDFILT2 Complex 2D median filter
if any(rem(d,2)~=1)
error('%s: Median filer window should have odd dimensions.',...
upper(mfilename));
end
dims = size(f);
if dims(1) < d(1)
error('%s: Median filer height is bigger than the image height.',...
upper(mfilename));
end
if dims(2) < d(2)
error('%s: Median filer width is bigger than the image width.',...
upper(mfilename));
end
fo = f;
whalf = floor(d(2)/2);
hhalf = floor(d(1)/2);
% Safely inside of the image
% The border values are taken from the input
for ii=(1+hhalf):(dims(1)-hhalf)
for jj=(1+whalf):(dims(2)-whalf)
neigh = sort(reshape(f(ii-hhalf:ii+hhalf,jj-whalf:jj+whalf),[],1));
fo(ii,jj) = neigh(ceil(end/2));
end
end
% Top rows
for ii=1:hhalf
for jj=(1+whalf):(dims(2)-whalf)
neigh = sort(reshape(f(1:ii+hhalf,jj-whalf:jj+whalf),[],1));
fo(ii,jj) = neigh(ceil(end/2));
end
end
% Top bottom rows
for ii=1:hhalf
for jj=(1+whalf):(dims(2)-whalf)
neigh = sort(reshape(f(end+1-(ii+hhalf):end,jj-whalf:jj+whalf),[],1));
fo(end+1-ii,jj) = neigh(ceil(end/2));
end
end
% The leftmost and rightmost collumns are not processed...
function [Ha,Hb]=parpeak(fc,Q,G,Fs)
% PARLSF Parametric Peaking filter
% Input parameters:
% fm : Cut-off frequency
% Q : Filter quality. Q=fc/B, where B is filter bandwidth.
% G : Gain in dB
% Fs : Sampling frequency
% Output parameters:
% Ha : Transfer function numerator coefficients.
% Hb : Transfer function denominator coefficients.
%
% For details see Table 5.3 in the Zolzer: Digital Audio Signal
% Processing, 2nd Edition, ISBN: 978-0-470-99785-7
Ha = zeros(3,1);
Hb = zeros(3,1);
%b0
Hb(1) = 1;
Ha(1) = 1;
K = tan(pi*fc/Fs);
if G>0
V0=10^(G/20);
den = 1 + K/Q + K*K;
% a0
Ha(1) = (1+V0*K/Q+K*K)/den;
% a1
Ha(2) = 2*(K*K-1)/den;
% a2
Ha(3) = (1-V0*K/Q+K*K)/den;
% b1
Hb(2) = 2*(K*K-1)/den;
% b2
Hb(3) = (1-K/Q+K*K)/den;
elseif G<0
V0=10^(-G/20);
den = 1 + V0*K/Q + K*K;
% a0
Ha(1) = (1+K/Q+K*K)/den;
% a1
Ha(2) = 2*(K*K-1)/den;
% a2
Ha(3) = (1-K/Q+K*K)/den;
% b1
Hb(2) = 2*(K*K-1)/den;
% b2
Hb(3) = (1-V0*K/Q+K*K)/den;
end