DSPRelated.com
   Code Snippets

DWT filter generator - formfilter.m

David Valencia ● October 28, 2010 ● Coded in Matlab 
% ------------------------------------------------------
% Title: formfilter.m
%
% Author: David Valencia
%
% Institution: UPIITA-IPN 2010
%
% for DSPrelated.com
% 
% Published in: http://www.dsprelated.com/showcode/12.php
%
% Description: This program receives 2 basic ortogonal
%               filters (one is high-pass and one is low-pass) 
%               the filtering level and the branch number.
%               As a result it returns the equivalent filter
%               of the specified branch.
% 
% Dependencies: upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revision: v1.0a
% - Commented and translated in English
%
% For more details on this code and its implementation
% see the following blog posts:
%
% http://www.dsprelated.com/showarticle/115.php
% http://www.dsprelated.com/showarticle/116.php
%
% ------------------------------------------------------

function [hx] = formfilter(n_stages,branch,h0,h1)
p = branch;

% Seed vector
hx = 0;
hx(1) = 1;

switch n_stages
    case 1
        % If it is a single stage filter
        % the branches are just the base filters
        if mod(branch,2) ~= 0
            hx = h0;
        else
            hx = h1;
        end
    case 2
        % For a 2 stage filter bank, the
        % equivalent filters are simply 
        % the convolution between the corresponding
        % base filters, and one of them is upsampled
        % The use of upsample2 is needed to avoid having
        % a large quantitiy of zeros at the end of the vector
        % which certainly difficults its usage to build a
        % convolution matrix.
        switch branch
            case 1
                hx = conv(h0,upsample2(h0,2));
            case 2
                hx = conv(h0,upsample2(h1,2));
            case 3
                hx = conv(h1,upsample2(h0,2));
            case 4
                hx = conv(h1,upsample2(h1,2));
            otherwise
                beep;
                fprintf('\nFor a 2 stage filter bank there can not be a fifth branch');
        end
        
    otherwise
        % For a n>2 stages filter bank, a more ellaborated
        % process must be made. A series of upsamplings and convolutions
        % are needed to get an equivalent vector.
        for i=0:n_stages-2
            q = floor(p /(2^(n_stages-1-i)));
            if (q == 1)
                hx = conv(hx,upsample2(h1,2^i));
            else
                hx = conv(hx,upsample2(h0,2^i));
            end
            p = mod(p,2^(n_stages-1-i));
        end
        
        % Depending on the parity of the branch number, the filter 
        % goes through a last convolution
        t = mod(branch,2);
        if(t == 1)
            hx = conv(hx,upsample2(h0,2^(n_stages-1)));
        else
            hx = conv(hx,upsample2(h1,2^(n_stages-1)));
        end
             
end

Discrete Wavelet Transform via linear convolution matrix

David Valencia ● October 28, 2010●2 comments ● Coded in Matlab 
% ----------------------------------------------------
% Title: Discrete Wavelet Transform
%        by means of a linear convolution matrix
%
% Author: David Valencia
% UPIITA IPN
%
% Posted at: http://www.dsprelated.com/showcode/11.php
%
% Description:
% Computes the Discrete Wavelet Transform 
% of n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% convolution matrix method (linear convolution)
%
% Dependencies: 
%
% formfilter.m
% http://www.dsprelated.com/showcode/12.php
% upsample2.m
% http://www.dsprelated.com/showcode/10.php
%
% Revisions:
%       v1.0a Commented and translated in English
%
% ----------------------------------------------------

close all; clear all; clc;

disp('== GENERALIZED MATRICIAL DWT ==')

%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db3';

% Obtain base filters
if(typeofbasis == 'b')
    [Rf,Df] = biorwavf(typbior);
    [h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
    [h0,h1,g0,g1] = wfilters(typor);
end;

%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.

N = 16; %Number of input samples
x = (1:N)';
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
dec_fact = n_branches; %Decimating factor

% L = Basic filters lenght (h0 รณ h1)
% N = Input vector length
% n_stages = number of stages, it generates (2^n_stages) branches

hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);

fprintf('\nLinear convolution matrix processing\n');

%% STEP 3: Build analysis matrices

% -------------------------------------------------------------
% High-pass analysis matrix (single scale) [linear convolution]
% As the DWT has a simple high-pass branch, the matrix is easy 
% to build with the base filter.

% Append zeros to each side of the vector to use it as a sliding
% part to build the matrix
movil_1 = [zeros(1,N-1),fliplr(h1),zeros(1,N-1)];
lm = length(movil_1);
n_rows = ceil((L+N-1)/2); %Compute the number of rows
W1 = zeros(n_rows,N); % Allocate memory for the high pass matrix
dec_factor = 2; %Decimate factor for the high-pass matrix

% Build matrix (with sequential shifting)
for i = 0:(n_rows-1)
    W1(i+1,:) = movil_1(lm-(i*dec_factor)-N+1 : lm-(i*dec_factor)); %Rama pasaaltas
end

disp('High-pass analysis matrix: ');
W1

disp('High-pass branch results');
y1 = W1 * x

% -------------------------------------------------------------
% Low-pass several stage analysis matrix [linear convolution]

% Get an equivalent filter to get its length
hx = formfilter(n_stages, 1, h0, h1);
Lhx = length(hx);

% Compute the number of rows per branch
rowsperfilter = ceil((N+Lhx)/dec_fact);

% Allocate memory for the low-pass matrix
W0 = zeros((n_branches/2)*rowsperfilter,N);

tic; %Start cronometer

% Build low pass filter
for i = 0:(n_branches/2)-1
    % Get the equivalent filter depending on the number of stages
    if n_stages < 3
        hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
    else
        hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
    end
    % Append zeros to the vector
    movil = [zeros(1,N-1),fliplr(hx),zeros(1,N-1)];
    lm = length(movil);
    % Shift and add the vector to the matrix
    for j = 0:rowsperfilter-1
        W0(i*rowsperfilter + j + 1,:) = movil(lm-(dec_fact*j)-N + 1 : lm-(dec_fact*j));
    end
end

disp('Low-pass many stage analysis matrix: ');
W0

disp('Low-pass filter bank results: ');
y0 = W0*x

%% STEP 4: Recover signal

disp('Recovered signal: ');
xn1 = (W1')*y1;
xn0 = (W0')*y0;
xn = xn1+xn0

%% STEP 5: Compute error
err = (x - xn).^2;
disp('Error :');
err = (sum(err)/N);
disp(err);

fprintf('\nProcessing time = %d seg \n',toc);
fprintf('\nProgram end >>> \n');

upsample2 - Upsample without extra overhead zeros

David Valencia ● October 27, 2010 ● Coded in Matlab 
% David Valencia
% Posted @ DSPRelated.com
% in: http://www.dsprelated.com/showcode/10.php

function [y] = upsample2(x,N)

lx = length(x);
y = 0;

% Put first sample in the first position of the output

y(1) = x(1);

for i=2:lx
    % Create a "chunk" (or piece) of vector, composed
	% with a series of zeros and the next sample from
	% the input
	chunk = [zeros(1,N-1),x(i)];
	
	% Append the chunk to the current output vector
    y = [y,chunk];
end

CIC Filter Analysis script

● October 27, 2010●6 comments ● Coded in Matlab 
%==========================================================================
% CIC Filter Analysis Script
%
% Compute CIC filter frequency response, compensation filter coefficients,
% and compensated frequency response for a given set of design paramenters.
%
%
%   Author:  Tom Chatt
%
%   Revision History:
%==========================================================================
%   09/23/2010              Tom Chatt                   Initial Version
%==========================================================================

clear all; close all; clc;

%=============================
% CIC Filter Design Parameters
%=============================
N = 3;                     % Sinc/Comb-Filter order
M = 1;                     % Stage delay (should be 1 or 2).
D = 64;                    % Decimation factor
fs = 10e6;                 % Pre-decimation sampling Frequency
fs_lo = fs / D;            % Post-decimation sampling frequency

%==============================
% Compensator Design Parameters
%==============================

B = 12; % Compensator coefficient bit width

% fir2 parameters
fc = fs/D/2;      % Ideal cutoff frequency
L = 16;           % Filter order - 1. Must be even. Final filter length is L + 1.
Fo = D*fc/fs;     % Normalized cutoff frequency

%=====================================================
% Generate and plot filters
%
% Floating point compensator filter coef. stored in h
% Fixed point compensator filter coef. stored in hz
%
%======================================================

% CIC Compensator using fir2
p = 2e3;                    % Granularity
s = 0.25 / p;               % Step size
fpass = 0:s:Fo;             % Pass band frequency samples
fstop = (Fo + s) : s : 0.5; % Stop band frequency samples
ff = [fpass fstop] * 2;     % Normalized frequency samples, 0 <= f <= 1

Mp = ones(1,length(fpass)); %% Pass band response; Mp(1)=1
Mp(2:end) = abs( M*D*sin(pi*fpass(2:end)/D)./sin(pi*M*fpass(2:end))).^N;
Mf = [Mp zeros(1,length(fstop))];
ff(end) = 1;
h = fir2(L,ff,Mf); %% Filter length L+1
h = h/max(h); %% Floating point coefficients
hz = round(h*power(2,B-1)-1)/(2^(B-1)); %% Fixed point coefficients

f = linspace(0,1,4*length(ff));
f_lo = [ff ff+1 ff+2 ff+3];
H = (D^-N*abs(D*M*sin(pi*M*D*f) ./ (pi*M*D*f)).^N);
HH = (D^-N*abs(D*M*sin(pi*M*ff) ./ (pi*M*ff)).^N);
HHH = (D^-N*abs(D*M*sin(pi*M*f_lo) ./ (pi*M*f_lo)).^N);

H_comp = abs(fft(hz,length(ff)));

%=========
% Plotting
%=========

figure('name','CIC Frequency Response','Numbertitle','off');
subplot(211);
plot(f*fs, db(H), 'b');
grid on;
title([{'Wideband frequency response of CIC filter'};{sprintf('D = %i, N = %i',D,N)}]);
xlabel('Frequency (Hz)');
ylabel('|H(e^{j\omega})| (dB)');
axis([0,fs,-1250,0]);

subplot(212);
plot(f_lo*fs_lo, db(HHH));
grid on;
title('Narrowband frequency response of CIC filter');
xlabel('Frequency (Hz)');
ylabel('|H(e^{j\omega})| (dB)');
axis([0,4*fs_lo,-200,5]);

figure('name','CIC with Compensator Filter','Numbertitle','off');
subplot(211);
hold on;
grid on;
plot(f_lo*fs_lo, db(HHH), 'b');
plot(f_lo*fs_lo, db([H_comp H_comp H_comp H_comp]), 'g');
plot(f_lo*fs_lo, db(HHH .* [H_comp H_comp H_comp H_comp]), 'r');
axis([0,4*fs_lo,-200,25])
title([{'CIC with Compensator FIR'};{sprintf('Compensator length = %i, Coefficient precision = %i',L,B)}]);
xlabel('Frequency (Hz)');
ylabel('|H(e^{j\omega})| (dB)');

subplot(212);
hold on;
grid on;
plot(ff*fs_lo, HH, 'b');
plot(ff*fs_lo, H_comp, 'g');
plot(ff*fs_lo, HH .* H_comp, 'r');
axis([0,1*fs_lo,0,1.1]);
title('Closeup of CIC + compensator passband');
xlabel('Frequency (Hz)');
ylabel('|H(e^{j\omega})|');
legend('CIC','Compensator Filter','Net Response','location','northeast');

Discrete Fourier Transform (DFT)

Miguel De Jesus Rosario ● October 27, 2010●11 comments ● Coded in Matlab 
function xdft=fdft(x)

N=length(x);

xdft(1,N)=0; % defining a 1xN sequence with all the elements equals to zero.

%*************THIS IS THE DFT IMPLEMENTATION**************
%For each value of k, the Summatory function is computed over the interval 0<n<N-1.
 
for k=0:N-1,   
  for n=0:N-1,
    xdft(k+1)=xdft(k+1) + x(n+1)*exp(-j*2*pi*k*n/N); 
  end
end

%NOTE: You may skip the use of for-loops by defining k and n as two N-points sequences. This is:

% k=0:N-1;  n=0:N-1; 

%The DFT is now computed in matrix form (as a matrix multiplication of the N-points Fourier Matrix and the discrete-time sequence x).

Linear and nonlinear phase low pass filter

● October 26, 2010●6 comments ● Coded in Matlab 
%% FIR & IIR filter
% Initializatoin

clc
clear all 
close all

%% FIR linear phase low pass filter

h=remez(60, [0 0.25 0.3 1], [1 1 0 0]);
fvtool(h,1)

%% IIR nonlinear phase low pass filter

[b a]=cheby2(20,35,0.3);
fvtool(b,a)

%% Filtering x(n)

x=zeros(1,200);
x(1:30)=1;

y1=filter(h,1,x);
y2=filter(b,a,x);

figure
plot(y1,'b','linewidth',2)
hold on
plot(y2,'r','linewidth',2)
xlabel('Samples')
ylabel('Amplitude')
title('Output of filters for input x(n)')
legend('FIR linear phase LPF output','IIR nonlinear phase LPF filter output')

Least-Squares silver bullets: Moore-Penrose Pseudoinverse (Example)

Markus Nentwig ● October 24, 2010●4 comments ● Coded in Matlab 
% *****************************************************************
% Using the Moore-Penrose Pseudoinverse to remove known disturbances from
% measured data
%
% all vectors are column vectors.
% *****************************************************************
close all; clear all;

% *****************************************************************
% Example: A measured signal has been polluted with power-line hum,
% at a fundamental frequency, 3rd and 5th harmonic, with unknown
% amplitude and phase. 
% *****************************************************************
m=1000; % number of measured points
fs=22000; % sampling rate
fHum=50; % known 50 Hz power line hum frequency
index=(1:m)';

% create power-line hum terms
hum=0.87654*sin(2*pi*index*fHum/fs+0.12345); 
hum3rd=0.65432*sin(2*pi*index*3*fHum/fs+0.45678); 
hum5th=0.3321*sin(2*pi*index*5*fHum/fs+0.78901); 

% create actual data
idealData=0.1*randn(m, 1);

% create measured data - deteriorated by hum
measuredData=idealData+hum+hum3rd+hum5th;

% *****************************************************************
% Processing to remove the hum
% *****************************************************************
% We know frequency and number of harmonics
% We DO NOT the amplitude and phase of each term.
% note: a*sin(ometaT)+b*cos(omegaT) is equivalent to c*sin(omegaT+d)
% where either [a, b] or [c, d] are the unknown parameters.

basis=[sin(2*pi*index*1*fHum/fs), cos(2*pi*index*1*fHum/fs), ...
       sin(2*pi*index*3*fHum/fs), cos(2*pi*index*3*fHum/fs), ...
       sin(2*pi*index*5*fHum/fs), cos(2*pi*index*5*fHum/fs)];

% *****************************************************************
% Moore-Penrose Pseudoinverse: Least-squares fit between the basis
% waveforms and the measured signal.
% *****************************************************************
leastSquaresCoeff=pinv(basis)*measuredData;
reconstructedHum=basis*leastSquaresCoeff;

correctedData=measuredData-reconstructedHum;

% *****************************************************************
% Plot measured signal and reconstructed hum
% *****************************************************************
figure(); hold on; 
plot(measuredData);
plot(reconstructedHum, 'r', 'LineWidth', 2);
legend({'measured data', 'reconstructed hum signal'});

% *****************************************************************
% Plot the remainder after subtracting the reconstructed hum.
% Compare with ideal result
% *****************************************************************
figure(); hold on;
plot(idealData,  'k', 'LineWidth', 2);
plot(correctedData);
plot(idealData-correctedData, 'r', 'LineWidth', 3);
legend({'ideal result', 'measured and corrected', 'error'})

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

Raised-cosine window for OFDM intersymbol transitions

Markus Nentwig ● October 21, 2010●6 comments ● Coded in Matlab 
% *************************************************************
% Raised cosine window for OFDM intersymbol transitions
% *************************************************************
% Use the window for the guard period between adjacent symbols
% to smoothen the waveform transition, limit sidelobes and prevent 
% unwanted emissions into adjacent frequency ranges
%
% reference:
% Data transmission by Frequency Division Multiplexing 
% using the Discrete Fourier Transform
% S.B. Weinstein and Paul M. Ebert
% IEEE transactions on Communications technology, No. 5,
% October 1971
% *************************************************************

% Length of the guard period in samples
n=64; 

index=linspace(0, 1, n+2);
index=index(2:end-1);

windowNextSymbol=(1-cos(pi*index))/2;
windowPreviousSymbol=fliplr(windowNextSymbol);
  
% plot

figure(); title('raised cosine window for OFDM symbol transitions');
xlabel('sample index into the guard period between symbols');
ylabel('window weight');
hold on; 
plot(windowPreviousSymbol, 'r+');
plot(windowNextSymbol, 'b+');
legend({'cyclic extension of symbol i-1', 'cyclic extension of symbol i'});

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'});