DSPRelated.com
   Code Snippets

Farrow resampler

Markus Nentwig ● August 12, 2011●4 comments ● Coded in Matlab 
% **************************************************************
% Vectorized Farrow resampler
% M. Nentwig, 2011
% Note: Uses cyclic signals (wraps around)
% **************************************************************
close all; clear all;

% inData contains input to the resampling process (instead of function arguments)
inData = struct();

% **************************************************************
% example coefficients. 
% Each column [c0; c1; c2; c3] describes a polynomial for one tap coefficent in fractional time ft [0, 1]:
% tapCoeff = c0 + c1 * ft + c2 * ft ^ 2 + c3 * ft ^ 3
% Each column corresponds to one tap. 
% the matrix size may be changed arbitrarily.
% 
% The example filter is based on a 6th order Chebyshev Laplace-domain prototype.
% **************************************************************
inData.cMatrix =[ -8.57738278e-3 7.82989032e-1 7.19303539e+000 6.90955718e+000 -2.62377450e+000 -6.85327127e-1 1.44681608e+000 -8.79147907e-1 7.82633997e-2 1.91318985e-1 -1.88573400e-1 6.91790782e-2 3.07723786e-3 -6.74800912e-3
 2.32448021e-1 2.52624309e+000 7.67543936e+000 -8.83951796e+000 -5.49838636e+000 6.07298348e+000 -2.16053205e+000 -7.59142947e-1 1.41269409e+000 -8.17735712e-1 1.98119464e-1 9.15904145e-2 -9.18092030e-2 2.74136108e-2
 -1.14183319e+000 6.86126458e+000 -6.86015957e+000 -6.35135894e+000 1.10745051e+001 -3.34847578e+000 -2.22405694e+000 3.14374725e+000 -1.68249886e+000 2.54083065e-1 3.22275037e-1 -3.04794927e-1 1.29393976e-1 -3.32026332e-2
 1.67363115e+000 -2.93090391e+000 -1.13549165e+000 5.65274939e+000 -3.60291782e+000 -6.20715544e-1 2.06619782e+000 -1.42159644e+000 3.75075865e-1 1.88433333e-1 -2.64135123e-1 1.47117661e-1 -4.71871047e-2 1.24921920e-2
] / 12.28;

% **************************************************************
% Create example signal
% **************************************************************
nIn = 50; % used for test signal generation only
if true
    % complex signal
    inData.signal = cos(2*pi*(0:(nIn-1)) / nIn) + cos(2*2*pi*(0:(nIn-1)) / nIn) + 1.5;
else
    % impulse response
    inData.signal = zeros(1, nIn); inData.signal(1) = 1; % 

    % inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 1) = 1; % enable to show constant c in first tap
    % inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(2, 1) = 1; % enable to show linear c in first tap
    % inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(3, 1) = 1; % enable to show quadratic c in first tap
    % inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 2) = 1; % enable to show constant c in second tap
end

% **************************************************************
% Resample to the following number of output samples
% must be integer, otherwise arbitrary
% **************************************************************
inData.nSamplesOut = floor(nIn * 6.28);

% **************************************************************
% Set up Farrow resampling
% **************************************************************
nSamplesIn = size(inData.signal, 2);
nSamplesOut = inData.nSamplesOut;
order = size(inData.cMatrix, 1) - 1; % polynomial order
nTaps = size(inData.cMatrix, 2); % number of input samples that contribute to one output sample (FIR size)

% pointer to the position in the input stream for each output sample (row vector, real numbers), starting at 0
inputIndex = (0:nSamplesOut-1) / nSamplesOut * nSamplesIn;

% split into integer part (0..nSamplesIn - 1) ...
inputIndexIntegerPart = floor(inputIndex);
% ... and fractional part [0, 1[
inputIndexFractionalPart = inputIndex - inputIndexIntegerPart;

% **************************************************************
% Calculate output stream
% First constant term (conventional FIR), then linear, quadratic, cubic, ... 
% **************************************************************
outStream = zeros(1, inData.nSamplesOut);
for ixOrder = 0 : order
    x = inputIndexFractionalPart .^ ixOrder;
    
    % **************************************************************
    % Add the contribution of each tap one-by-one
    % **************************************************************
    for ixTap = 0 : nTaps - 1
        c = inData.cMatrix(ixOrder+1, ixTap+1);
        % index of input sample that contributes to output via the current tap
        % higher tap index => longer delay => older input sample => smaller data index
        dataIx = inputIndexIntegerPart - ixTap;
        % wrap around
        dataIx = mod(dataIx, nSamplesIn);
        % array indexing starts at 1
        dataIx = dataIx + 1;
        delayed = inData.signal(dataIx);
        % for each individual output sample (index in row vector), 
        % - evaluate f = c(order, tapindex) * fracPart .^ order 
        % - scale the delayed input sample with f
        % - accumulate the contribution of all taps
        % this implementation performs the operation for all output samples in parallel (row vector)
        outStream = outStream + c * delayed .* x;
    end % for ixTap
end % for ixOrder

% **************************************************************
% plot
% **************************************************************
xIn = linspace(0, 1, nSamplesIn + 1); xIn = xIn(1:end-1);
xOut = linspace(0, 1, nSamplesOut + 1); xOut = xOut(1:end-1);

figure(); grid on; hold on;
stem(xIn, inData.signal, 'k+-');
plot(xOut, outStream, 'b+-');
legend('input', 'output');
title('Farrow resampling. Signals are cyclic');

Continuous-time equivalent of a discrete-time impulse response

Markus Nentwig ● August 12, 2011 ● Coded in Matlab 
% ************************************************************
% * Construct a continuous-time impulse response based on a discrete-time filter
% ************************************************************
close all; clear all;

% ************************************************************
% * example filter
% * f = [0 0.7 0.8 1]; a = [1 1 0 0];
% * b = remez(45, f, a);
% * h = b .';
% Illustrates rather clearly the difficulty of "interpolating" (in a geometric sense, 
% via polynomials, splines etc) between impulse response samples
% ************************************************************    
h = [2.64186706e-003 2.50796828e-003 -4.25450455e-003 4.82982358e-003 -2.51252769e-003 -2.52706568e-003 7.73569146e-003 -9.30398382e-003 4.65100223e-003 5.17152459e-003 -1.49409856e-002 1.76001904e-002 -8.65545521e-003 -9.78478603e-003 2.82103205e-002 -3.36173643e-002 1.68597129e-002 2.01914744e-002 -6.17486493e-002 8.13362871e-002 -4.80981494e-002 -8.05143565e-002 5.87677665e-001 5.87677665e-001 -8.05143565e-002 -4.80981494e-002 8.13362871e-002 -6.17486493e-002 2.01914744e-002 1.68597129e-002 -3.36173643e-002 2.82103205e-002 -9.78478603e-003 -8.65545521e-003 1.76001904e-002 -1.49409856e-002 5.17152459e-003 4.65100223e-003 -9.30398382e-003 7.73569146e-003 -2.52706568e-003 -2.51252769e-003 4.82982358e-003 -4.25450455e-003 2.50796828e-003 2.64186706e-003];

% ************************************************************
% zero-pad the impulse response of the FIR filter to 5x its original length
% (time domain)
% The filter remains functionally identical, since appending zero taps changes nothing
% ************************************************************
timePaddingFactor = 5;
n1 = timePaddingFactor * size(h, 2);
nh = size(h, 2);
nPad = n1 - nh;
nPad1 = floor(nPad / 2);
nPad2 = nPad - nPad1;
h = [zeros(1, nPad1), h, zeros(1, nPad2)];
hOrig = h;

% ************************************************************
% Determine equivalent Fourier series coefficients (frequency domain)
% assume that the impulse response is bandlimited (time-discrete signal by definition) 
% and periodic (period length see above, can be extended arbitrarily)
% ************************************************************
h = fft(h); % Fourier transform time domain to frequency domain

% ************************************************************
% Evaluate the Fourier series on an arbitrarily dense grid 
% this allows to resample the impulse response on an arbitrarily dense grid
% ************************************************************
% zero-pad Fourier transform
% ideal band-limited oversampling of the impulse response to n2 samples
n2 = 10 * n1;

h = [h(1:ceil(n1/2)), zeros(1, n2-n1), h(ceil(n1/2)+1:end)];
h = ifft(h); % back to time domain
% Note: One may write out the definition of ifft() and evaluate the exp(...) term at an 
% arbitrary time to acquire a true continuous-time equation.

% numerical inaccuracy in (i)fft causes some imaginary part ~1e-15
assert(max(abs(imag(h))) / max(abs(real(h))) < 1e-14);
h = real(h);

% scale to maintain amplitude level
h = h * n2 / n1;

% construct x-axis [0, 1[
xOrig = linspace(-timePaddingFactor/2, timePaddingFactor/2, size(hOrig, 2) + 1); xOrig = xOrig(1:end-1);
x = linspace(-timePaddingFactor/2, timePaddingFactor/2, size(h, 2) + 1); x = x(1:end-1);

% ************************************************************
% Plot 
% ************************************************************
% ... on a linear scale

% find points where original impulse response is defined (for illustration)
ixOrig = find(abs(hOrig) > 0);

figure(); grid on; hold on;
stem(xOrig(ixOrig), hOrig(ixOrig), 'k');
plot(x, h, 'b'); 
legend('discrete-time FIR filter', 'continuous-time equivalent');
title('equivalent discrete- and continuous-time filters');
xlabel('time (relative to discrete time IR duration)'); ylabel('h(t)');

% ... and again on a logarithmic scale
myEps = 1e-15;
figure(); grid on; hold on;
u = plot(xOrig(ixOrig), 20*log10(abs(hOrig(ixOrig)) + myEps), 'k+'); set(u, 'lineWidth', 3);
plot(x, 20*log10(abs(h) + myEps), 'b');
legend('discrete-time FIR filter', 'continuous-time equivalent');
title('equivalent discrete- and continuous-time filters');
xlabel('time (relative to discrete time IR duration)'); ylabel('|h(t)| / dB');

First order RC Filter in the Digital Domain

Jeff T ● August 7, 2011●1 comment ● Coded in Matlab 
function [b,a] = rc_filter(R, C, Fs, filter_type)
% Returns equivalent IIR coefficients for an analog RC filter
%
% Usage:     [B,A] = RC_FILTER(r, c, fs, type);
%
%             R is the resistance value (in ohms)
%             C is the capacitance value (in farrads)
%             FS is the digital sample rate (in Hz)
%             type is a character string defining filter type
%                 Choices are: 'high' or 'low'
%
% Highpass filters have the analog structure:
%
%
%                  | | 
%     Vi  o--------| |----------+---------o  Vo
%                  | |          |
%                     C         | 
%                              --- 
%                              | |  R
%                              | | 
%                              ---   
%                               |
%                               |
%         o---------------------+---------o
%                              GND
%
%
% Lowpass filters have the analog structure:
%
%
%                  |-----| 
%     Vi  o--------|     |------+---------o  Vo
%                  |-----|      |
%                     R         | 
%                             ----- C
%                             ----- 
%                               | 
%                               |
%         o---------------------+---------o
%                              GND
%
% This function uses a pre-calculated equation for both of these circuits
% that only requires the resistance and capacitance value to get a true
% digital filter equivalent to a basic analog filter.
%
% The math behind these equations is based off the basic bilinear transform
% technique that can be found in many DSP textbooks.  The reference paper
% for this function was "Conversion of Analog to Digital Transfer 
% Functions" by C. Sidney Burrus, page 6.
%
% This is also the equivalent of a 1st order butterworth with a cuttoff
% frequency of Fc = 1/(2*pi*R*C);
%
% Author: sparafucile17 07/01/02
%

% Verify that cutoff of the analog filter is below Nyquist
if(  (1/(2*pi*R*C))  >  (Fs/2)  )
    error('This analog filter cannot be realized with this sample rate');
end

% Default to allpass if invalid type is selected
b = [ 1 0 ];
a = [ 1 0 ];

% Constants
RC = R * C;
T  = 1 / Fs;

% Analog Cutoff Fc
w = 1 / (RC);

% Prewarped coefficient for Bilinear transform
A = 1 / (tan((w*T) / 2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
%            s
% T(s) =  -------
%          s + 1
%
%
% using Bilinear transform of
%
%             1          ( 1 - z^-1 )
% s -->  -----------  *  ------------
%         tan(w*T/2)     ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(strcmp(filter_type,'high'))
    
    b(1) = (A)     / (1 + A);
    b(2) = -b(1);
    a(2) = (1 - A) / (1 + A);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
%            1
% T(s) =  -------
%          s + 1
%
%
% using Bilinear transform of
%
%             1          ( 1 - z^-1 )
% s -->  -----------  *  ------------
%         tan(w*T/2)     ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(strcmp(filter_type,'low'))
    
    b(1) = (1)     / (1 + A);
    b(2) = b(1);
    a(2) = (1 - A) / (1 + A);

end

Exponential Audio Unmute

Jeff T ● August 7, 2011 ● Coded in Matlab 
function [y] = signal_unmute(x, index, duration, Fs)
% Return the input vector with a unmute that occurs at a specific
% index and with an exponential ramp-up to reduce "pop" sounds.
% The output will start off at -100dB gain and end at 0dB gain.
%
% Usage: y = SIGNAL_UNMUTE(x, index, duration, Fs);
%
%        X is your one-dimensional input array
%        INDEX is where in the input signal you want the unmute to begin
%        DURATION is how long (in seconds) to exponentially ramp up the
%           input signal. 100ms is recommended.
%        FS is the sample rate
%
% Example: 
%        You want to unmute your signal at 0.5sec with a duration of 100ms 
%        and with a 48k Fs. (signal is unmuted completely at 0.6sec)
%        y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003

% Input must have some length
if(length(x) == 1)
    error('ERROR: input signal must have more than one element');
end

% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
    error('ERROR: Input must be one-dimensional');
end

% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
    error(['There are not enough samples in X to complete the unmute.  '  ...
           'Either change the mute duration or move the index back.' ]);
end

% Flip vector (temporarily)
if((size(x, 2) ~= 1))
    x = x';
    flip_back = true;
else
    flip_back = false;
end

% Calculate exponential coefficient
dB_atten  = -100; %What do we consider "muted" to be in decibels
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));

% Build the Gain array
gain = [zeros(index, 1); ones((length(x)-index), 1);];
b = [ coeff  0];
a = [ 1 (-1*(1-coeff))];
gain  = filter(b,a,gain);

% Apply Mute (gain) to the input signal
y = gain .* x;

% Flip the vector (if required)
if(flip_back == true);
    y = y';
end

Exponential Audio Mute

Jeff T ● August 7, 2011 ● Coded in Matlab 
function [y] = signal_mute(x, index, duration, Fs)
% Return the input vector with a mute that occurs at a specific
% index and with an exponential ramp-down to reduce "pop" sounds.
% The output will start off at 0dB gain and end at -100dB gain.
%
% Usage: y = SIGNAL_MUTE(x, index, duration, Fs);
%
%        X is your one-dimensional input array
%        INDEX is where in the input signal you want the mute to begin
%        DURATION is how long (in seconds) to exponentially ramp down the
%           input signal. 100ms is recommended.
%        FS is the sample rate
%
% Example: 
%        You want to mute your signal at 0.5sec with a duration of 100ms 
%        and with a 48k Fs. (mute complete at 0.6sec)
%        y = signal_mute(x, 24000, 0.1, 48000);
%
% Author: sparafucile17 7/29/2003

% Input must have some length
if(length(x) == 1)
    error('ERROR: input signal must have more than one element');
end

% This function only supports one-dimensional arrays
if((size(x, 2) ~= 1) && (size(x, 1) ~= 1))
    error('ERROR: Input must be one-dimensional');
end

% Make sure there are enough samples to complete the mute
if(length(x) < (index + duration*Fs))
    error(['There are not enough samples in X to complete the mute.  '  ...
           'Either change the mute duration or move the index back.' ]);
end

% Flip vector (temporarily)
if((size(x, 2) ~= 1))
    x = x';
    flip_back = true;
else
    flip_back = false;
end

% Calculate exponential coefficient
dB_atten  = -100; %How much attenuation will be at time: index + duration
decayrate = -dB_atten / duration;
coeff = 1.0 - 10^(-decayrate/(20.0*Fs));

% Build the Gain array
gain = [ones(index, 1); zeros((length(x)-index), 1);];
b = [ coeff  0];
a = [ 1 (-1*(1-coeff))];
gain  = filter(b,a,gain, [1]);

% Remove minor overshoot at the beginning of gain vector
gain(find(gain > 1)) = 1;

% Apply Mute (gain) to the input signal
y = gain .* x;

% Flip the vector (if required)
if(flip_back == true);
    y = y';
end

Image denoising: Threshold calculation for vishushrink method

Senthilkumar ● July 30, 2011 ● Coded in Matlab 
%function to calculate the threshold using Visushrink denoising method
function T = Visu_threshold(X)
    [m,n] = size(X);
    M = m*n;
    T = sqrt(2*log(M));

Image denoising: using vishushrink method

Senthilkumar ● July 30, 2011 ● Coded in Matlab 
function [soft_X1,SOFT_PSNR] = Vishu_soft(X,Y)
%function used to denoise a noisy image using vishuShrink method
%One -level decomposition
[CA,CH,CV,CD] = dwt2(Y,'haar');
%Call the function to calculate the threshold
T1 = Visu_threshold(CD);
% Call the function to perfom soft shrinkage
de_CH = soft(CH,T1);
de_CV = soft(CV,T1);
de_CD = soft(CD,T1);
%
%Two -level decomposition
[CA1,CH1,CV1,CD1] = dwt2(CA,'haar');
%Call the function to calculate the threshold
T2 = Visu_threshold(CD1);
% Call the function to perfom soft shrinkage
de_CH1 = soft(CH1,T2);
de_CV1 = soft(CV1,T2);
de_CD1 = soft(CD1,T2);
% % CA1 = soft1(CA1,T2);
%
%
%Reconstruction for soft shrinkage
X2 = idwt2(CA1,de_CH1,de_CV1,de_CD1,'haar');
X1 = idwt2(X2,de_CH,de_CV,de_CD,'haar');

SOFT_PSNR = PSNR(X,X1);
soft_X1 = uint8(X1);

Image denoise: Threshold for sureshrink method

Senthilkumar ● July 30, 2011●1 comment ● Coded in Matlab 
function thr = sureshrink(CD,T)
%function used to calculate threshold using sureshrink method
CD = CD(:)';
n = length(CD);
sx2 = sort(abs(CD)-T).^2;  % sorting in descending order
b  = cumsum(sx2);           %cumulative sum
risks = (n-(2*(1:n))+b)/n;
[risk,best] = min(risks);
thr = sqrt(sx2(best));

Image denoise: denoising using soft threshold

Senthilkumar ● July 30, 2011 ● Coded in Matlab 
function y = soft(x,T)
%function used to denoise a noisy image with given soft threshold
%x = noisy image
%T = threshold
y = max(abs(x) - T, 0);
y = y./(y+T) .* x;

Image denoising: Peak signal to Noise ratio calculation

Senthilkumar ● July 30, 2011●1 comment ● Coded in Matlab 
function A = PSNR(G,H)
%G = original image
%H =  denoised image
  error = G - H;
  decibels = 20*log10((255*255)/(sqrt(mean(mean(error.^2)))));
  disp(sprintf('PSNR = %5.2f db',decibels)) 
  A = decibels;