DSP Code Snippets

Miguel De Jesus Rosario ● December 14, 2010●1 comment ● Coded in Matlab

% VectorGoertzel    Goertzel's Algorithm filter bank.
%   Realization of the Goertzel's Algorithm to compute the Nonuniform DFT 
%   of a signal(a column vector named signalw) of length Nw with sampling 
%   frecuency fs at the desired frecuencies contained in vector f. The 
%   function returns the NDFT magnitude in a vector with the same length of f.

function xk=Gfilterbank(signalw,f,fs,Nw)

% Inititialization of the different variables 
n=2;
signalw=signalw(1:Nw);

cost=cos(2*pi*f/fs);
sint=sin(2*pi*f/fs);
L=length(cost);

y1(1:L)=0;
y2(1:L)=0;
signalw=[0 0 signalw]; %Signal is delayed by two samples

% Goertzel Feedback Algorithm

while((n-2) < Nw)
    n=n+1;
    xnew(1:L)=signalw(n);
    y=xnew+2*cost.*y1-y2;
    y2=y1;
    y1=y;
end

% Goertzel Forward Algorithm

rey=y1-y2.*cost;
imy=y2.*sint;

% Magnitude Calculation

xk=abs(rey+imy*j)';

Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

Julius Orion Smith III ● November 8, 2010●7 comments ● Coded in Matlab

NZ = 1;      % number of ZEROS in the filter to be designed
NP = 4;      % number of POLES in the filter to be designed
NG = 10;     % number of gain measurements
fmin = 100;  % lowest measurement frequency (Hz)
fmax = 3000; % highest measurement frequency (Hz)
fs = 10000;  % discrete-time sampling rate
Nfft = 512;  % FFT size to use 
df = (fmax/fmin)^(1/(NG-1)); % uniform log-freq spacing
f = fmin * df .^ (0:NG-1);   % measurement frequency axis

% Gain measurements (synthetic example = triangular amp response):
Gdb = 10*[1:NG/2,NG/2:-1:1]/(NG/2); % between 0 and 10 dB gain

% Must decide on a dc value.
% Either use what is known to be true or pick something "maximally
% smooth".  Here we do a simple linear extrapolation:
dc_amp = Gdb(1) - f(1)*(Gdb(2)-Gdb(1))/(f(2)-f(1));

% Must also decide on a value at half the sampling rate.
% Use either a realistic estimate or something "maximally smooth".
% Here we do a simple linear extrapolation. While zeroing it
% is appealing, we do not want any zeros on the unit circle here.
Gdb_last_slope = (Gdb(NG) - Gdb(NG-1)) / (f(NG) - f(NG-1));
nyq_amp = Gdb(NG) + Gdb_last_slope * (fs/2 - f(NG));

Gdbe = [dc_amp, Gdb, nyq_amp]; 
fe = [0,f,fs/2];
NGe = NG+2;

% Resample to a uniform frequency grid, as required by ifft.
% We do this by fitting cubic splines evaluated on the fft grid:
Gdbei = spline(fe,Gdbe); % say `help spline'
fk = fs*[0:Nfft/2]/Nfft; % fft frequency grid (nonneg freqs)
Gdbfk = ppval(Gdbei,fk); % Uniformly resampled amp-resp

figure(1);
semilogx(fk(2:end-1),Gdbfk(2:end-1),'-k'); grid('on'); 
axis([fmin/2 fmax*2 -3 11]);
hold('on'); semilogx(f,Gdb,'ok');
xlabel('Frequency (Hz)');   ylabel('Magnitude (dB)');
title(['Measured and Extrapolated/Interpolated/Resampled ',...
       'Amplitude Response']);

Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies

S = 10 .^ (Sdb/20); % convert to linear magnitude
s = ifft(S); % desired impulse response
s = real(s); % any imaginary part is quantization noise
tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
disp(sprintf(['Time-limitedness check: Outer 20%% of impulse ' ...
              'response is %0.2f %% of total rms'],tlerr));
% = 0.02 percent
if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
  error('Increase Nfft and/or smooth Sdb');
end

figure(2);
plot(s, '-k'); grid('on');   title('Impulse Response');
xlabel('Time (samples)');   ylabel('Amplitude');

c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
% Check aliasing of cepstrum (in theory there is always some):
caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
disp(sprintf(['Cepstral time-aliasing check: Outer 20%% of ' ...
    'cepstrum holds %0.2f %% of total rms'],caliaserr));
% = 0.09 percent
if caliaserr>1.0 % arbitrary limit
  error('Increase Nfft and/or smooth Sdb to shorten cepstrum');
end
% Fold cepstrum to reflect non-min-phase zeros inside unit circle:
% If complex:
% cf = [c(1), c(2:Ns-1)+conj(c(Nfft:-1:Ns+1)), c(Ns), zeros(1,Nfft-Ns)];
cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
Cf = fft(cf); % = dB_magnitude + j * minimum_phase
Smp = 10 .^ (Cf/20); % minimum-phase spectrum

Smpp = Smp(1:Ns); % nonnegative-frequency portion
wt = 1 ./ (fk+1); % typical weight fn for audio
wk = 2*pi*fk/fs;
[B,A] = invfreqz(Smpp,wk,NZ,NP,wt);
Hh = freqz(B,A,Ns);

figure(3);
plot(fk,db([Smpp(:),Hh(:)])); grid('on');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Frequency Response');
legend('Desired','Filter');

NCO bit true model in a single equation

kaz - ● October 23, 2011 ● Coded in Matlab

%NCO bit true model
clear all; close all;

%example NCO parameters
n = 20000;                     %number of test samples
A = 2^15 - 1;                  %max amplitude, amplitude resolution 16 bit signed
M = 2^32;                      %NCO accumulator phase resolution
inc = round(M*.001);           %NCO accumulator phase increment 
k = 2^12;                      %lut phase resolution (lut size), may equal M

%single equation for NCO bit true model
y1 = round(A*sin(2*pi*round(k*mod((0:n-1)*inc/M,M))/k));      

%%%%%%%%%%%%%%%%%%%%%%%%% equivalent code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lsb = log2(M) - log2(k);                 %LSBs discarded when addressing
lut = round(A*sin(2*pi*(0:k-1)/k));      %lut, one cycle sine data

ptr = 0; 
addr = 0;
for i = 1:n 
    y2(i) = lut(addr+1);          %add 1 for matlab LUT
    ptr = mod(ptr + inc, M);      %phase accumulator(ramp)
    addr = round(ptr/2^lsb);      %discard LSBs
    addr(addr >= k) = addr - k;   %check address overflow
end

%compare
plot(y1);hold;
plot(y2,'r--')
plot(y1 - y2,'g--')

Complex filter with three multipliers per tap

kaz - ● October 22, 2011 ● Coded in Matlab

%random input
I_in = randn(1,1024);
Q_in = randn(1,1024);

%random test filter
h = randn(1,20)+j*randn(1,20);    

%split up into 3 subfilters
h1 = real(h)+imag(h);         %extra adder
h2 = real(h);
h3 = real(h)-imag(h);         %extra subtractor

%apply subfilter
ch1 = filter(h1,1,I_in);
ch2 = filter(h2,1,Q_in-I_in); %extra subtractor
ch3 = filter(h3,1,Q_in);

%combine outputs
yr = round((ch3-ch2)*2^15);
yi = round((ch1+ch2)*2^15);
IQ_out = complex(yr,yi);

%direct computation
ref = round(2^15*filter(h,1,complex(I_in,Q_in)));
plot(IQ_out - ref);

Beampattern of a Linear Array of Antennas (Array Processing)

Miguel De Jesus Rosario ● November 30, 2010●1 comment ● Coded in Matlab

% The function NDFT2 computes the sinc pattern of a linear array of 
% sensors. The function receives three input: the interelement spacing d,
% the array weigths a, and the zero padding desired Npad.

function NDFT2(d, a, Npad)

k  = -Npad/2 : Npad/2-1; % index
u  = pi*k/Npad; % u is a vector with elements in the range of -pi/2 to pi/2
uk = sin(u);
n  = 0:Npad-1; % n is an index
m  = 1:Npad; % m is an index

% Wavenumber K = 2*pi/landa 
% d = is a fraction or a multiple of lambda.
% v = K*d*uk(m).
  
v = 2*pi*d*uk(m)';
Dn(m,n+1) = exp(j*v*n);

puk = Dn * a'; % Computing the Beampattern

% Plot of the Beampatterns
figure(1); subplot(2,1,1); plot(uk,20*log10(abs(puk)));grid on; xlabel('sin(theta)'); ylabel('Magnitude (dB)')
axis([-1 1 -100 0])
subplot(2,1,2); plot(uk,puk);grid on; xlabel('sin(theta)'); ylabel('Magnitude')
axis([-1 1 -1 1])
warning off;
% Plot the beampattern as a polar graph
figure(2); polar(u',puk); hold on; polar(-u',puk);hold off         

%  Function pattern()
%
%     The function pattern() computes and plots the beampattern of a 
%     Linear Array of sensors. This function has three inputs: the number of elements in
%     the array, the pointing direction of the beam pattern (an angular
%     direction in radians), and a constant spacing between the elements of 
%     the array (as fraction of the wavelenght(lambda)of the transmitted  
%     signal.  The optimum spacing between elements of the array is 0.5. 
%     You must also select any of the windows presented in the menu. 
%     Windowing techniques are used to reduced the sidelobes of the pattern.   
%     The list of available windows is the following:
%
%           @bartlett       - Bartlett window.
%           @barthannwin    - Modified Bartlett-Hanning window. 
%           @blackman       - Blackman window.
%           @blackmanharris - Minimum 4-term Blackman-Harris window.
%           @bohmanwin      - Bohman window.
%           @chebwin        - Chebyshev window.
%           @flattopwin     - Flat Top window.
%           @gausswin       - Gaussian window.
%           @hamming        - Hamming window.
%           @hann           - Hann window.
%           @kaiser         - Kaiser window.
%           @nuttallwin     - Nuttall defined minimum 4-term Blackman-Harris window.
%           @parzenwin      - Parzen (de la Valle-Poussin) window.
%           @rectwin        - Rectangular window.
%           @tukeywin       - Tukey window.
%           @triang         - Triangular window.
%
%     For example:      pattern(21, pi/4, 0.5);
%           Plots the beampattern of an linear array with 21 elements equally spaced  
%           at a half of the wavelenght(lambda/2), and a pointing
%           direction of 45 degrees. For uniform arrays use a rectangular 
%           window (rectwin).

function pattern(array_number, angular_direction, spacing)

close all
clc

N = array_number;
delta = angular_direction;
d = spacing;

Npad=1024;
n=0:N-1;

delta = 2*pi*d*sin(delta);
an=1/N*exp(j*n*delta);

for i=0:500000,
    
option = menu('Choose the desired Window', 'Bartlett', 'Barthannwin', 'Blackman', 'Blackmanharris', 'Bohmanwin', 'Chebwin', 'Flattopwin', 'Gausswin', 'Hamming', 'Hann', 'Kaiser', 'Nuttallwin', 'Parzenwin', 'Rectwin', 'Tukeywin', 'Triang', 'Exit'); 

switch option
    
    case 1
        close all
        clear a;
        a=an;
        a = a.*bartlett(N)';
        a(Npad)=0;
      
        NDFT2(d, a, Npad);
        
        
    case 2
        close all
        clear a;
        a=an;
        a = a.*barthannwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
        
    case 3
        close all 
        clear a;
        a=an;
        a = a.*blackman(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 4
        close all  
        clear a;
        a=an;
        a = a.*blackmanharris(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 5
        close all 
        clear a;
        a=an;
        a = a.*bohmanwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 6
        close all 
        clear a;
        a=an;
        a = a.*chebwin(N,40)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 7
        close all
        clear a;
        a=an;
        a = a.*flattopwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);
        
    case 8
        close all 
        clear a;
        a=an;
        a = a.*gausswin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 9
        close all  
        clear a;
        a=an;
        a = a.*hamming(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 10
        close all 
        clear a;
        a=an;
        a = a.*hann(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 11
         close all
         clear a;
        a=an;
        a = a.*kaiser(N,1)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 12
         close all
         clear a;
        a=an;
        a = a.*nuttallwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 13
         close all 
         clear a;
        a=an;
        a = a.*parzenwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 14
         close all 
         clear a;
        a=an;
        a = a.*rectwin(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 15
         close all  
         clear a;
        a=an;
        a = a.*tukeywin(N,0)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 16
        close all
        clear a;
        a=an;
        a = a.*triang(N)';
        a(Npad)=0;

        NDFT2(d, a, Npad);

    case 17
        
        break;
              
    end
    
end

Adding a Controlled Amount of Noise to a Noise-Free Signal

Rick Lyons ● January 3, 2012●21 comments ● Coded in Matlab

function [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
%   SNR_Set(x, Desired_SNR_dB) returns the real-valued 
%   input 'Signal' contaminated with normally-distributed, 
%   zero-mean, random noise.  The signal-to-noise ratio   
%   (SNR in dB) of the output 'Noisy_Signal' signal is  
%   controlled by the input 'Desired_SNR_dB' argument measured 
%   in dB.

%   Example:
% 	Npts = 128;                  % Number of time samples
% 	n = 0:Npts-1;                % Time-domain index
% 	Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% 	Desired_SNR_dB = 3;      % Set SNR of output 'Noisy_Signal' to +3 dB
% 	[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
%   Author: Richard Lyons [December 2011]
%******************************************

Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one

Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
	%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));

K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10);  % Scale factor

New_Noise = sqrt(K)*Noise; % Change Noise level
	%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
	%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))

Noisy_Signal = Signal + New_Noise;

'SNR_Set()' Function Test Code:
%  Filename SNR_Set_test.m
%
%  Tests the 'SNR_Set()" function.  Adds a predefined 
%  amount of random noise to a noise-free signal such that  
%  the noisy signal has a desired signal-to-noise ratio (SNR).
%
%  Author: Richard Lyons [December 2011]

clear, clc

% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5;  % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal

Desired_SNR_dB = 3 % Set desired SNR in dB

[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal

% Plot original and 'noisy' signals
    figure(1)
    subplot(2,1,1)
    plot(n, Signal, '-bo', 'markersize', 2)
    title('Original Signal')
    grid on, zoom on
    subplot(2,1,2)
    plot(n, Noisy_Signal, '-bo', 'markersize', 2)
    title('Noisy Signal')
    xlabel('Time-samples')
    grid on, zoom on

% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
    figure(2)
    plot(Spec_Mag, '-bo', 'markersize', 2)
    title('Spec Mag of Noisy Signal')
    xlabel('Freq-samples'), ylabel('Linear')
    grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)

Echo Filter

Miguel De Jesus Rosario ● December 1, 2010●3 comments ● Coded in Matlab

% function echo_filter:
%     1) computes and returns the impulse response h of the LTI system (echo filter)
%     2) Filter the input signal.  The output signal is yecho. 
%
% INPUTS:
% signal = input signal
% times = vector containing the time delays (in seconds) of the repetitions(echoes)
% attenuations = vector containing the respective attenuations of the echoes
% fs = sampling frequency of the input signal.
%
% OUTPUTS:
% yecho = output signal (passed through the echo filter)
% h = impulse response of the echo filter.

function [yecho, h] = echo_filter(signal, times, attenuations, fs)

h(1)=1; % This is the first coefficient of the echo filter

% Computing the impulse response h
for i=1:length(times),
    samples(i) = times(i)*fs; % Calculating the sample-times (N = t*fs)
    h(floor(samples(i))) = attenuations(i); % Impulse response coeficients
end

% #########################################################################
% You may use the following implementation instead of the illustrated
% above:

%       samples = times*fs;
%       h(floor(samples)) = attenuations;

% In this case, the implementation is done without a for-loop.
% #########################################################################

yecho = filter(h,1,signal(:,1)); % filtering of the input signal results in a signal with echoes

% You may test the filter using the Matlab signal "splat" (write load splat
% in the Matlab command window).  Use function sound() in order to
% listening both, the input and the output, signals.

Generate ideal QPSK..64QAM symbol error rates (reference result)

Markus Nentwig ● October 21, 2010 ● Coded in Matlab

% ******************************************************
% Symbol error rate of Quadrature Amplitude Modulation
% ******************************************************
% uncoded, additive white Gaussian noise channel
% Implements [1] Figure 5.2-16 on a wider SNR range
% [1] John G. Proakis
%     "Digital Communications"
%      4th edition, McGraw-Hill", 2001
% http://www.amazon.com/Digital-Communications-John-Proakis/dp/0072321113

% Note that different definitions are in common use:
% a) {symbol, bit} errors
% b) signal-to-noise ratio per {symbol, bit}

close all; clear all;
SNR_perSymbol_dB=[-3:0.1:42];

% definition: Below [1] 5.2-78
SNR_perSymbol=10.^(SNR_perSymbol_dB/10);

plotPerBitX=[]; plotY=[];

% number of constellation points (M-ary QAM)
for M=[4, 16, 64, 256, 1024]

  % the energy of each symbol is used to transmit log2(M) bits
  SNR_perBit=SNR_perSymbol/log2(M);
   
  % [1] 5.2-78 argument of Q(...)
  Qarg=sqrt(3/(M-1)*SNR_perSymbol);
 
  % [1] 2.1-98
  Q=1/2*erfc(Qarg/sqrt(2));
   
  % [1] eq. 5.2-77
  % probability of error for PAM per quadrature signal
  PsqrtM=2*(1-1/sqrt(M))*Q;
 
  % [1] eq. 5.2-79
  PM=1-(1-PsqrtM).^2;

  plotPerBitX=[plotPerBitX; 10*log10(SNR_perBit)];
  plotY=[plotY; PM];
end

figure(1);
h = semilogy(SNR_perSymbol_dB', plotY'); grid on;
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per symbol (dB)'); ylabel('symbol error rate');
xlim([0, 40]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});

figure(2);
h = semilogy(plotPerBitX', plotY'); grid on; ylabel('symbol error rate');
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per bit (dB)');
xlim([-3, 30]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});

Shelving Filter Design

Jeff T ● May 29, 2011●1 comment ● Coded in Matlab

function [b, a]  = shelving(G, fc, fs, Q, type)

%
% Derive coefficients for a shelving filter with a given amplitude and
% cutoff frequency.  All coefficients are calculated as described in 
% Zolzer's DAFX book (p. 50 -55).  
%
% Usage:     [B,A] = shelving(G, Fc, Fs, Q, type);
%
%            G is the logrithmic gain (in dB)
%            FC is the center frequency
%            Fs is the sampling rate
%            Q adjusts the slope be replacing the sqrt(2) term
%            type is a character string defining filter type
%                 Choices are: 'Base_Shelf' or 'Treble_Shelf'
%
% Author:    sparafucile17 08/22/05
%

%Error Check
if((strcmp(type,'Base_Shelf') ~= 1) && (strcmp(type,'Treble_Shelf') ~= 1))
    error(['Unsupported Filter Type: ' type]);
end

K = tan((pi * fc)/fs);
V0 = 10^(G/20);
root2 = 1/Q; %sqrt(2)

%Invert gain if a cut
if(V0 < 1)
    V0 = 1/V0;
end

%%%%%%%%%%%%%%%%%%%%
%    BASE BOOST
%%%%%%%%%%%%%%%%%%%%
if(( G > 0 ) & (strcmp(type,'Base_Shelf')))
   
    b0 = (1 + sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
    b1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*K + K^2);
    b2 = (1 - sqrt(V0)*root2*K + V0*K^2) / (1 + root2*K + K^2);
    a1 =                (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
    a2 =             (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%    BASE CUT
%%%%%%%%%%%%%%%%%%%%
elseif (( G < 0 ) & (strcmp(type,'Base_Shelf')))
    
    b0 =             (1 + root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
    b1 =                (2 * (K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
    b2 =             (1 - root2*K + K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);
    a1 =             (2 * (V0*K^2 - 1) ) / (1 + root2*sqrt(V0)*K + V0*K^2);
    a2 = (1 - root2*sqrt(V0)*K + V0*K^2) / (1 + root2*sqrt(V0)*K + V0*K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE BOOST
%%%%%%%%%%%%%%%%%%%%
elseif (( G > 0 ) & (strcmp(type,'Treble_Shelf')))

    b0 = (V0 + root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
    b1 =             (2 * (K^2 - V0) ) / (1 + root2*K + K^2);
    b2 = (V0 - root2*sqrt(V0)*K + K^2) / (1 + root2*K + K^2);
    a1 =              (2 * (K^2 - 1) ) / (1 + root2*K + K^2);
    a2 =           (1 - root2*K + K^2) / (1 + root2*K + K^2);

%%%%%%%%%%%%%%%%%%%%
%   TREBLE CUT
%%%%%%%%%%%%%%%%%%%%

elseif (( G < 0 ) & (strcmp(type,'Treble_Shelf')))

    b0 =               (1 + root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
    b1 =                  (2 * (K^2 - 1) ) / (V0 + root2*sqrt(V0)*K + K^2);
    b2 =               (1 - root2*K + K^2) / (V0 + root2*sqrt(V0)*K + K^2);
    a1 =             (2 * ((K^2)/V0 - 1) ) / (1 + root2/sqrt(V0)*K + (K^2)/V0);
    a2 = (1 - root2/sqrt(V0)*K + (K^2)/V0) / (1 + root2/sqrt(V0)*K + (K^2)/V0);

%%%%%%%%%%%%%%%%%%%%
%   All-Pass
%%%%%%%%%%%%%%%%%%%%
else
    b0 = V0;
    b1 = 0;
    b2 = 0;
    a1 = 0;
    a2 = 0;
end

%return values
a = [  1, a1, a2];
b = [ b0, b1, b2];

Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Markus Nentwig ● October 24, 2010●12 comments ● Coded in Matlab

% *****************************************************************
% Audio interpolation using Cubic Hermite Spline (Catmull-Rom)
%
% The audio file samples describe a continuous-waveform. 
% The task is to interpolate between samples, for example in sample rate
% conversion or variable-speed playback.
% The cubic spline approach is one particular trade-off between 
% quality and efficiency.
% 
% For an efficient C-implementation, refer to fluid_rvoice_dsp.c
% in the fluidsynth project on sourceforge.
% Based on work by Olli Niemitalo:
% http://yehar.com/blog/wp-content/uploads/2009/08/deip.pdf page 11

% ******************************************************************
% build coefficient matrix
% ******************************************************************
% ideally, the interpolation coefficients would be recalculated for any possible
% value of ptrFrac below.
% example input file: http://www.elisanet.fi/mnentwig/webroot/common/in.wav
% Note: This example is written for clarity, not speed.
nPoly=128;
x=((0:nPoly-1)/nPoly)'; % column vector

c1=(x .* (-0.5 + x .* (1 - 0.5 .* x)));
c2=(1.0 + x .* x .* (1.5 .* x - 2.5));
c3=(x .* (0.5 + x .* (2.0 - 1.5 .* x)));
c4=(0.5 .* x .* x .* (x - 1.0));

% ******************************************************************
% load input data
% ******************************************************************
[y, fs, bits]=wavread('in.wav');
[nSamples, nChan]=size(y);
inputdata=transpose(y(:, 1)); % first channel if stereo

% prepare output data storage
outIndex=1; outputData=[]; nAlloc=0;

% ******************************************************************
% "ptr" indicates a position in the input data. 
% Loop through the whole range until end-of-data.
% ******************************************************************
ptr=2;
while ptr < nSamples-2

    % "rate" determines the number of input samples per output sample
    % rate=1.1; % faster constant-speed playback
    % variable playback speed (arbitrary example)
    rate=1+0.3*sin(outIndex/40000);;

    % ******************************************************************
    % Interpolate samples
    % ******************************************************************
    % ptr is the index into the input data, needed for the next output sample.
    % Usually, it lies betwen points, for example ptr=1234.56

    % integer part (1234)
    ptrBase=floor(ptr);
    % fractional part [0..1[, for example 0.56
    ptrFrac=ptr-ptrBase;

    % look up the filter for the fractional part
    fIndex=floor(ptrFrac*nPoly)+1;

    % input samples
    s1=inputdata(ptrBase-1);
    s2=inputdata(ptrBase);
    s3=inputdata(ptrBase+1);
    s4=inputdata(ptrBase+2);

    % performance optimization. Reallocate some extra space.
    if outIndex>=nAlloc
        nAlloc=nAlloc+2000;
        outputData(nAlloc)=0;
    end

    % ******************************************************************
    % Calculate output sample.
    % ******************************************************************
    outputData(outIndex)=s1*c1(fIndex)+s2*c2(fIndex)+s3*c3(fIndex)+s4*c4(fIndex);

    % that's it. Next output sample
    outIndex=outIndex+1;

    % advance the index in the source data
    ptr=ptr+rate;
end

% performance optimization part II, chop extra memory
outputData=outputData(1:outIndex-1)';

% ******************************************************************
% Write output to file
% ******************************************************************
% wavwrite('out.wav', outputData, fs, bits); % OCTAVE version
wavwrite(outputData, fs, bits, 'out.wav'); % MATLAB version

Previous 8 91011 12 13 Next

Image Quantizer Demo

Ron ● April 1, 2011 ● Coded in Matlab

clf;
for (truncate = 0:8)
x = imread('aqua.jpg');
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
x_t = (x/2^truncate) * 2^truncate;
%if above line doesn't work your MATLAB is old! use next line:
% x_t = uint8(double(uint8((double(x)/2^truncate))) * 2^truncate);
imshow(x_t);
subplot(2,2,3);
imshow(uint8(abs(double(x_t) - double(x))));
pause(0.5);
end

Signal-to-Distortion SDR Ratio

Ron ● April 1, 2011 ● Coded in Matlab

[f, fs, nbits] = wavread('input.wav');
[f2, fs2, nbits2] = wavread('input-truncated.wav');

top = 1/size(f,1)*(sum(f.^2))
bottom = 1/size(f2,1)*(sum(f2.^2))

Y = 10*log10(top/bottom)

RGB (JPEG) Image Separator

Ron ● April 1, 2011 ● Coded in Matlab

clf;
x = imread('input.jpg');
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);

subplot(2,2,1);
imshow(x);
title('Original');
subplot(2,2,2);
imshow(red);
title('Red');
subplot(2,2,3);
imshow(green);
title('Green');
subplot(2,2,4);
imshow(blue);
title('Blue');

Image Quantizer

Ron ● April 1, 2011 ● Coded in Matlab

%number of bits to truncate
tred = 4;
tgreen = 2;
tblue = 3;

clf;
x = imread('input.jpg');
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);

red_t = (red/2^tred) * 2^tred;
green_t = (green/2^tgreen) * 2^tgreen;
blue_t = (blue/2^tblue) * 2^tblue;

x_t(:,:,1) = red_t;
x_t(:,:,2) = green_t;
x_t(:,:,3) = blue_t;

imshow(x);
pause;
imshow(x_t)

size = (24 - tred - tgreen - tblue) / 24

2D Quantizer with Error Output

Ron ● April 1, 2011 ● Coded in Matlab

clf;
truncate = 5;
x = imread('input.jpg');
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
x_t = (x/2^truncate) * 2^truncate;
%if above line doesn't work your MATLAB is old! use next line:
% x_t = uint8(double(uint8((double(x)/2^truncate))) * 2^truncate);
imshow(x_t);
subplot(2,2,3);
imshow(uint8(abs(double(x_t) - double(x))));

WAV Quantizer with Error Output

Ron ● April 1, 2011 ● Coded in Matlab

truncate = 6;
[f, fs, nbits] = wavread('in.wav');
f_int = int16(f*2^(nbits-1-truncate));
f_back_to_float = (double(f_int))/2^(nbits-1-truncate);
f_diff = f - f_back_to_float;
wavwrite(f_back_to_float, fs, nbits, 'out.wav');
wavwrite(f_diff, fs, nbits, 'error.wav');

Quantization Demo

Ron ● April 1, 2011 ● Coded in Matlab

a = 10;
x = 0:2*pi/1000:2*pi;
y = a*sin(x);
y_int = int8(y);

subplot(3,1,1);
plot(x,y);
grid;
axis([0 2*pi -a*1.2 a*1.2]);
title('Continuous');

subplot(3,1,2);
plot(x,y_int);
grid;
axis([0 2*pi -a*1.2 a*1.2]);
title('Quantized');

subplot(3,1,3);
plot(x,abs(double(y_int)-y));
grid;
axis([0 2*pi 0 1.2/2]);
title('Error');

WAV Subsample

Ron ● April 1, 2011 ● Coded in Matlab

factor = 8;
[X, fs, nbits] = wavread('in.wav');
Y = downsample(X,factor);
wavwrite (Y, fs/factor, nbits, 'out.wav');

2D Image Subsample

Ron ● April 1, 2011 ● Coded in Matlab

A = imread('in.jpg');

M =   [ 1 4  6  4  1;
        4 16 24 16 4;
        6 24 36 24 6;
        4 16 24 16 4;
        1 4  6  4  1 ] / 256;

B = uint8(conv2(double(A),M));
imwrite(B,'out.jpg','JPG');

2D Image FIR Subsample

Ron ● April 1, 2011 ● Coded in Matlab

A = imread('in.jpg');

M =   [ 1 4  6  4  1;
        4 16 24 16 4;
        6 24 36 24 6;
        4 16 24 16 4;
        1 4  6  4  1 ] / 256;

Y = filter2(M,A);
B = uint8(conv2(double(Y),M));

imwrite(B,'out.jpg','JPG');

Previous 8 91011 12 13 Next

Goertzel Filterbank to the Implementation of a Nonuniform DFT

Fitting Filters to Measured Amplitude Response Data Using invfreqz in Matlab

NCO bit true model in a single equation

Complex filter with three multipliers per tap

Beampattern of a Linear Array of Antennas (Array Processing)

Adding a Controlled Amount of Noise to a Noise-Free Signal

Echo Filter

Generate ideal QPSK..64QAM symbol error rates (reference result)

Shelving Filter Design

Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Image Quantizer Demo

Signal-to-Distortion SDR Ratio

RGB (JPEG) Image Separator

Image Quantizer

2D Quantizer with Error Output

WAV Quantizer with Error Output

Quantization Demo

WAV Subsample

2D Image Subsample

2D Image FIR Subsample

Sign in

Search code

Blogs - Hall of Fame

Free PDF Downloads

Quick Links

About DSPRelated.com

Social Networks

The Related Media Group