DSP Code Snippets

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

Fixed Point FIR Filter

Shawn Stevenson ● December 3, 2010●1 comment ● Coded in C

#include <stdio.h>
#include <stdint.h>

//////////////////////////////////////////////////////////////
//  Filter Code Definitions
//////////////////////////////////////////////////////////////

// maximum number of inputs that can be handled
// in one function call
#define MAX_INPUT_LEN   80
// maximum length of filter than can be handled
#define MAX_FLT_LEN     63
// buffer to hold all of the input samples
#define BUFFER_LEN      (MAX_FLT_LEN - 1 + MAX_INPUT_LEN)

// array to hold input samples
static int16_t insamp[ BUFFER_LEN ];

// FIR init
void firFixedInit( void )
{
    memset( insamp, 0, sizeof( insamp ) );
}

// store new input samples
int16_t *firStoreNewSamples( int16_t *inp, int length )
{
    // put the new samples at the high end of the buffer
    memcpy( &insamp[MAX_FLT_LEN - 1], inp,
            length * sizeof(int16_t) );
    // return the location at which to apply the filtering
    return &insamp[MAX_FLT_LEN - 1];
}

// move processed samples
void firMoveProcSamples( int length )
{
    // shift input samples back in time for next time
    memmove( &insamp[0], &insamp[length],
            (MAX_FLT_LEN - 1) * sizeof(int16_t) );
}

// the FIR filter function
void firFixed( int16_t *coeffs, int16_t *input, int16_t *output,
       int length, int filterLength )
{
    int32_t acc;     // accumulator for MACs
    int16_t *coeffp; // pointer to coefficients
    int16_t *inputp; // pointer to input samples
    int n;
    int k;

    // apply the filter to each input sample
    for ( n = 0; n < length; n++ ) {
        // calculate output n
        coeffp = coeffs;
        inputp = &input[n];
        // load rounding constant
        acc = 1 << 14;
        // perform the multiply-accumulate
        for ( k = 0; k < filterLength; k++ ) {
            acc += (int32_t)(*coeffp++) * (int32_t)(*inputp--);
        }
        // saturate the result
        if ( acc > 0x3fffffff ) {
            acc = 0x3fffffff;
        } else if ( acc < -0x40000000 ) {
            acc = -0x40000000;
        }
        // convert from Q30 to Q15
        output[n] = (int16_t)(acc >> 15);
    }
}

//////////////////////////////////////////////////////////////
//  Test program
//////////////////////////////////////////////////////////////

// bandpass filter centred around 1000 Hz
// sampling rate = 8000 Hz
// gain at 1000 Hz is about 1.13

#define FILTER_LEN  63
int16_t coeffs[ FILTER_LEN ] =
{
 -1468, 1058,   594,   287,    186,  284,   485,   613,
   495,   90,  -435,  -762,   -615,   21,   821,  1269,
   982,    9, -1132, -1721,  -1296,    1,  1445,  2136,
  1570,    0, -1666, -2413,  -1735,   -2,  1770,  2512,
  1770,   -2, -1735, -2413,  -1666,    0,  1570,  2136,
  1445,    1, -1296, -1721,  -1132,    9,   982,  1269,
   821,   21,  -615,  -762,   -435,   90,   495,   613,
   485,  284,   186,   287,    594, 1058, -1468
};

// Moving average (lowpass) filter of length 8
// There is a null in the spectrum at 1000 Hz

#define FILTER_LEN_MA   8
int16_t coeffsMa[ FILTER_LEN_MA ] =
{
    32768/8, 32768/8, 32768/8, 32768/8,
    32768/8, 32768/8, 32768/8, 32768/8
};

// number of samples to read per loop
#define SAMPLES   80

int main( void )
{
    int size;
    int16_t input[SAMPLES];
    int16_t output[SAMPLES];
    int16_t *inp;
    FILE   *in_fid;
    FILE   *out_fid;
    FILE   *out_fid2;

    // open the input waveform file
    in_fid = fopen( "input.pcm", "rb" );
    if ( in_fid == 0 ) {
        printf("couldn't open input.pcm");
        return;
    }

    // open the output waveform files
    out_fid = fopen( "outputFixed.pcm", "wb" );
    if ( out_fid == 0 ) {
        printf("couldn't open outputFixed.pcm");
        return;
    }
    out_fid2 = fopen( "outputFixedMa.pcm", "wb" );
    if ( out_fid == 0 ) {
        printf("couldn't open outputFixedMa.pcm");
        return;
    }

    // initialize the filter
    firFixedInit();

    // process all of the samples
    do {
        // read samples from file
        size = fread( input, sizeof(int16_t), SAMPLES, in_fid );
        // store new samples in working array
        inp = firStoreNewSamples( input, size );

        // apply each filter
        firFixed( coeffs, inp, output, size, FILTER_LEN );
        fwrite( output, sizeof(int16_t), size, out_fid );

        firFixed( coeffsMa, inp, output, size, FILTER_LEN_MA );
        fwrite( output, sizeof(int16_t), size, out_fid2 );

        // move processed samples
        firMoveProcSamples( size );
    } while ( size != 0 );

    fclose( in_fid );
    fclose( out_fid );
    fclose( out_fid2 );

    return 0;
}

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

Template to develop audio DSP using the (OSS) API + GNU/Linux

Gabriel Rivas ● February 6, 2011 ● Coded in C

#include <stdio.h>
#include <fcntl.h>
#include <sys/soundcard.h>

#define BUF_SIZE 2

/*buffers for sound device*/
unsigned char devbuf[BUF_SIZE];

/*File descriptor for device*/
int audio_fd;

/*The sound card is known as the DSP*/
char DEVICE_NAME[]="/dev/dsp";

/*Device Settings*/
int format;
int channels;
int fs;
int len;
int frag;
int devcaps;

int main()
{
    unsigned int temp;
    unsigned int yout;

    /*Samples format is 16 bit unsigned*/
    format = AFMT_U16_LE;
    /*1 Channel, MONO*/
    channels = 1;
    /*Sampling Rate is 16 KHz*/
    fs = 16000;
    /*This is a parameter used to set the DSP for low latencies*/
    frag = 0x00020004;

	/******************************************************
          Set parameters in sound card device
    *****************************************************/

	/*Open sound card device*/
	if ((audio_fd = open(DEVICE_NAME,O_RDWR, 0)) == -1) {
		/* Open of device failed */
		perror(DEVICE_NAME);
		exit(1);
	}	

	if (ioctl (audio_fd, SNDCTL_DSP_SETFRAGMENT, &frag) == -1)
	{
	    perror ("SNDCTL_DSP_SETFRAGMENT");
	    exit (-1);
	}

	/*Set audio format*/
	if (ioctl(audio_fd, SNDCTL_DSP_SETFMT, &format) == -1) {
		/* fatal error */
		perror("SNDCTL_DSP_SETFMT");
		exit(1);
	}
	/*Set number of channels*/
	if (ioctl(audio_fd, SNDCTL_DSP_CHANNELS, &channels) == -1) {
		/* Fatal error */
		perror("SNDCTL_DSP_CHANNELS");
		exit(1);
	}
	
	/*Set sampling rate*/
	if (ioctl(audio_fd, SNDCTL_DSP_SPEED, &fs)==-1) {
		/* Fatal error */
		perror("SNDCTL_DSP_SPEED");
		exit(1);
	}

	if (ioctl(audio_fd, SNDCTL_DSP_SYNC) == -1) {
		/* fatal error */
		perror("SNDCTL_DSP_SYNC");
		exit(1);
	}	

	/******************************************************
          This is the infinite loop to:
          1. Read a sample
          2. Process the sample
          3. Write the sample to soundcard output
    *****************************************************/
    while(1) {
        /* This is a blocking read function, the application will wait
           until the sound card captures a sample
        */
        if ((len = read(audio_fd,devbuf, sizeof(devbuf))) == -1) {
	        perror("audio read");
	        exit(1);
	    }

        /* We can pass this variable to a temporary value, 
           in this case as unsigned 16 value, and then use this value 
           for processing
        */  
	    temp = (devbuf[0])|(devbuf[1]<<8);

        /* In this case no processing is done so the value is passed to the output. You can also use this sample to build an audio file, or send it trought a communications interface, etc.
        */
	    yout = temp;

	    devbuf[0] = yout&0xFF;
	    devbuf[1] = (yout>>8)&0xFF;
		
        /* This is the way the output samples are sent to the soundcard            
        */
        if ((len = write(audio_fd,devbuf,sizeof(devbuf)) == -1)) {
	        perror("audio write");
	        exit(1);
        }
    }	
    return 0;
}

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

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

OFDM symbol generator

kaz - ● April 4, 2011●3 comments ● Coded in Matlab

%bandwidth = nCar/nFFT *Fs

clear all;
nData = 200000;                %vector length
ModulationType = '64QAM';
nFFT = 2048;                   %fft size
nCar = 1200;                   %number of active carriers
nPrefix = 144;
Fs = 10;                       %sampling rate
       
%compute number of ofdm symbols  
nSymbols = ceil(nData/(nFFT+nPrefix));
   
%generate random data symbols
switch ModulationType                              
    case 'QPSK',  alphabet = [-1,1];
    case '16QAM', alphabet = [-3,-1,1,3];
    case '64QAM', alphabet = [-7,-5,-3,-1,1,3,5,7];
end
   
tx = [];
for k = 1:nSymbols    
   QAM = complex(randsrc(nCar,1,alphabet),randsrc(nCar,1,alphabet));
 
   fft_in = zeros(1,nFFT);
   fft_in(1:nCar/2) = QAM(1:nCar/2);
   fft_in(end-nCar/2+1:end) = QAM(nCar/2+1:end);
   fft_out = ifft(fft_in);

   
   %concatenate ofdm symbol plus its prefix into tx vector
   prefix = fft_out(nFFT-nPrefix+1:nFFT);          
   tx = [tx prefix fft_out];
end

[P, F] = pwelch(tx, hann(2^16), 0, 2^16, Fs);
P = fftshift(P);
F = F - max(F)/2;
plot(F,10*log10(P));

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

Fast audio convolution

Hernán Ordiales ● November 19, 2010 ● Coded in Python

#!/usr/bin/python
# -*- coding: utf-8 -*-

#
# Copyright (C) 2006 by Hernán Ordiales
# <audiocode@uint8.com.ar>
#
# 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 2 of the License, or
# (at your option) any later version.
#

# Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
# Use: fast_conv.py source.wav impulse_response.wav output.wav

from wav_array import *
from FFT import fft, inverse_fft
from pylab import size
import sys

def nextpow2( L ):
	N = 2
	while N < L: N = N * 2
	return N
	
def fast_conv_vect( x, h ):
	# searches for the amount of points required to perform the FFT
	L = size(h) + size(x) - 1 # linear convolution length
	N = nextpow2( L )
	# Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
	# common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))

	# FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
	H = fft( h, N ) # Fourier transform of the impulse
	X = fft( x, N ) # Fourier transform of the input signal
	
	Y = H * X # spectral multiplication
	y = inverse_fft( Y ) # time domain again
	return y

source = sys.argv[1]
impulse = sys.argv[2]
conv_output = sys.argv[3]

clip_factor = 1.01 # clip factor default value

[ h1, Fs1, h_bits ] = wavread( impulse ) # impulse response
[ x1, Fs2, x_bits ] = wavread( source ) # file to process

if Fs1 == Fs2 : # if sample rates are the same
	print "Processing..."
	y1 = fast_conv_vect( x1, h1 ).real # takes the real part to avoid a too small complex part (around e-18)

	# audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
	y1 = y1/( max(y1)*clip_factor ) # to avoid clipping
	
	wavwrite( y1, Fs1, conv_output )

	print "Output file:", conv_output
	print "Convolution success."
else:
	print "Error: files has different sample rate."

Random Number Generator

● February 27, 2011●2 comments ● Coded in C for the SHARC

// Random number generator that sounds pleasant in audio algorithms.
// The argument and return value is a 30-bit (positive nonzero) integer.
// This is Duttaro's 30-bit pseudorandom number generator, p.124 J.AES, vol 50 no 3 March 2002; 
// I found this 30-bit pseudorandom number generator sounds far superior to the 31-bit 
// pseudorandom number generator presented in the same article.
// To make result a signed fixed-point value, do not shift left into the sign bit;
// instead, subtract 0x20000000 then multiply by a scale factor.
#define NextRand(rr) (((((rr >> 16) ^ (rr >> 15) ^ (rr >> 1) ^ rr) & 1) << 29) | (rr >> 1))

Previous 18 192021 22 23 Next

Analog IIR filter to Digital IIR Filter - Backward Difference

Senthilkumar ● March 21, 2011●3 comments ● Coded in Scilab

//Program To convert analog filter into digital filter using
//mapping = (z-(z^-1))/T - Backward Difference  
clear all;
clc;
close;
s = poly(0,'s');
H = 1/((s+0.1)^2+9)
T =1;//Sampling period T = 1 Second
z = poly(0,'z');
Hz = horner(H,(1/T)*(z-(z^-1)))

Frequency response of window functions for FIR Filter

Senthilkumar ● March 21, 2011●2 comments ● Coded in Scilab

//Program to find and plot the frequency response of
//(1) Hanning window (2)Hamming window for M = 11 
clear all;
close;
clc
M = 11;
for n = 1:M
  h_hann_11(n) = 0.5-0.5*cos(2*%pi*(n-1)/(M-1));
  h_hamm_11(n) = 0.54-0.46*cos(2*%pi*(n-1)/(M-1));
end
subplot(2,1,1)
[h_hann_11_M,fr]=frmag(h_hann_11,512);
h_hann_11_M = 20*log10(h_hann_11_M./max(h_hann_11_M));
plot2d(fr,h_hann_11_M,2);
title('Frequency Response 0f Hanning window Filter length M =11')
subplot(2,1,2)
[h_hamm_11_M,fr]=frmag(h_hamm_11,512);
h_hamm_11_M = 20*log10(h_hamm_11_M./max(h_hamm_11_M));
plot2d(fr,h_hamm_11_M,2);
title('Frequency Response of Hamming window Filter length M =11')

Window Functions for FIR Filter design

Senthilkumar ● March 21, 2011 ● Coded in Scilab

//Program to generate different window functions
//For FIR Filter design based on windowing techniques
clear all;
close;
clc
M =11 ;
for n = 1:M
  h_Rect(n) = 1;
  h_hann(n) = 0.5-0.5*cos(2*%pi*(n-1)/(M-1));
  h_hamm(n) = 0.54-0.46*cos(2*%pi*(n-1)/(M-1));
  h_balckmann(n) = 0.42-0.5*cos(2*%pi*n/(M-1))+0.08*cos(4*%pi*n/(M-1));
end
plot2d(1:M,[h_Rect,h_hann,h_hamm,h_balckmann],[2,5,7,9]); 
legend(['Rectangular Window';'Hanning';'Hamming';'Balckmann']);
title('Window Functions for Length M = 11')

Interpolation or upsampling in scilab

Senthilkumar ● March 21, 2011 ● Coded in Scilab

//Multirate Signal Processing in scilab
//Upsampling a sinusoidal signal by a factor of 2
clear;
clc;
n = 0:%pi/200:2*%pi;
x = sin(%pi*n);   //original signal
upsampling_x = zeros(1,2*length(x)); //upsampled by a factor of 2
upsampling_x(1:2:2*length(x)) = x;
subplot(2,1,1)
plot(1:length(x),x);
xtitle('original singal')
subplot(2,1,2)
plot(1:length(upsampling_x),upsampling_x);
xtitle('Upsampled Signal by a factor of 2');

Decimation or Downsampling in scilab

Senthilkumar ● March 21, 2011●2 comments ● Coded in Scilab

//Multirate Signal Processing in scilab
//Downsampling a sinusoidal signal by a factor of 2
clear;
clc;
n = 0:%pi/200:2*%pi;
x = sin(%pi*n);   //original signal
downsampling_x = x(1:2:length(x)); //downsampled by a factor of 2
subplot(2,1,1)
plot(1:length(x),x);
xtitle('original singal')
subplot(2,1,2)
plot(1:length(downsampling_x),downsampling_x);
xtitle('Downsampled Signal by a factor of 2');

Discrete Plot in scilab

Senthilkumar ● March 21, 2011●1 comment ● Coded in Scilab

//This Program Illustrates the discrete plot in scilab
//using plot2d3 function
clear;
clc;
close;
a =1.5; 
n =1:10;  
x = (a)^n;
a=gca();
a.thickness = 2;
plot2d3('gnn',n,x) 
xtitle('Graphical Representation of Exponentially Increasing Signal','n','x[n]');

Continuous Time Fourier Transform of a Square Waveform

Senthilkumar ● March 21, 2011 ● Coded in Scilab

// Continuous Time Fourier Transform and Frequency Response of a Square Waveform
// x(t)= A, from -T1 to T1
clear;
clc;
close;
// CTS Signal
A =1;    //Amplitude
Dt = 0.005;
T1 = 4;  //Time in seconds
t = -T1:Dt:T1;
for i = 1:length(t)
  xt(i) = A;
end
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1;        //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
xt = xt';
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = real(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
//
subplot(2,1,1);
a = gca();
a.data_bounds=[-4,0;4,2];
a.y_location ="origin";
plot2d2(t,xt);
xlabel('t in msec.');
title('Continuous Time Signal x(t)')
subplot(2,1,2);
a = gca();
a.y_location ="origin";
a.x_location ="origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds');
title('Continuous-time Fourier Transform  X(jW)')

Continuous Time Fourier Transform of a Exponential Signal

Senthilkumar ● March 21, 2011 ● Coded in Scilab

//Continuous Time Fourier Transform of a
//Continuous Time Exponential Signal x(t)= exp(-A*t)u(t), A>0
//And Plotting its Magnitude response and phase response
clear;
clc;
close;
// Analog Signal
A =1;    //Amplitude
Dt = 0.005;
t = 0:Dt:10;
xt = exp(-A*t);
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1;        //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = abs(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
[XW_Phase,db] = phasemag(XW);
XW_Phase = [-mtlb_fliplr(XW_Phase),XW_Phase(2:1001)];
//Plotting Continuous Time Signal
figure
a = gca();
a.y_location = "origin";
plot(t,xt);
xlabel('t in sec.');
ylabel('x(t)')
title('Continuous Time Signal')
figure
//Plotting Magnitude Response of CTS
subplot(2,1,1);
a = gca();
a.y_location = "origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds---> W');
ylabel('abs(X(jW))')
title('Magnitude Response (CTFT)')
//Plotting Phase Reponse of CTS
subplot(2,1,2);
a = gca();
a.y_location = "origin";
a.x_location = "origin";
plot(W,XW_Phase*%pi/180);
xlabel('                         Frequency in Radians/Seconds---> W');
ylabel('                                                <X(jW)')
title('Phase Response(CTFT) in Radians')

CTFS coefficients of a periodic square waveform

Senthilkumar ● March 21, 2011 ● Coded in Scilab

//CTFS coefficients of a periodic square waveform 
//x(t) = 1, |t|<T1, and 0, T1<|t|<T/2
clear;
close;
clc;
T =4;
T1 = T/4;
t = -T1:T1/100:T1;
Wo = 2*%pi/T;
xt =ones(1,length(t)); //Square Wave Generation
//Finding the Fourier Series Coefficients
for k =0:5
  C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
  a(k+1) = xt*C(k+1,:)'/length(t);
  if(abs(a(k+1))<=0.1) 
    a(k+1)=0;
  end
end
a =a';
a_conj = real(a(:))-sqrt(-1)*imag(a(:));
ak = [a_conj($:-1:1)',a(2:$)];
k = 0:5;
k = [-k($:-1:1),k(2:$)];
Spectrum_ak = (1/2)*real(ak);
//Plotting the Continuous Time Signal
figure
a = gca();
a.y_location = "origin";
a.x_location = "origin";
a.data_bounds=[-2,0;2,2];
plot2d2(t,xt,5)
poly1 = a.children(1).children(1);
poly1.thickness = 3; 
title('x(t)')
xlabel('                                                       t')
//Plotting the Spectral coefficients of CT Signal
figure
a = gca();
a.y_location = "origin";
a.x_location = "origin";
plot2d3('gnn',k,Spectrum_ak,5)
poly1 = a.children(1).children(1);
poly1.thickness = 3; 
title('abs(ak)')
xlabel('                                                       k')
//Result
//Spectrum_ak  =
//0.0633127  0. -0.1055559 0. 0.3167197 0.5 0.3167197 0. -0.1055559    0.    0.0633127

Continuous Time Fourier Series of sine signal

Senthilkumar ● March 21, 2011 ● Coded in Scilab

//Continuous Time Fourier Series Spectral Coefficients of 
//a periodic sine signal x(t) = sin(Wot)
clear;
close;
clc;
t = 0:0.01:1;
T = 1;
Wo = 2*%pi/T;
xt = sin(Wo*t);
for k =0:5
  C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
  a(k+1) = xt*C(k+1,:)'/length(t); //fourier series is done
  if(abs(a(k+1))<=0.01) 
    a(k+1)=0;
  end
end
a =a';
ak = [-a($:-1:1),a(2:$)];
disp(ak,'Continuous Time Fourier Series Coefficients are:')
//Result
//Continuous Time Fourier Series Coefficients are:   
// column  1 to 11
// 0 0 0 0 -1.010D-18+0.4950495i 0 1.010D-18-0.4950495i 0 0 0 0

Previous 18 192021 22 23 Next

Audio interpolation (Catmull-Rom/Cubic Hermite Spline)

Fixed Point FIR Filter

CIC Filter Analysis script

Template to develop audio DSP using the (OSS) API + GNU/Linux

Shelving Filter Design

First order RC Filter in the Digital Domain

OFDM symbol generator

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

Fast audio convolution

Random Number Generator

Analog IIR filter to Digital IIR Filter - Backward Difference

Frequency response of window functions for FIR Filter

Window Functions for FIR Filter design

Interpolation or upsampling in scilab

Decimation or Downsampling in scilab

Discrete Plot in scilab

Continuous Time Fourier Transform of a Square Waveform

Continuous Time Fourier Transform of a Exponential Signal

CTFS coefficients of a periodic square waveform

Continuous Time Fourier Series of sine signal

Sign in

Search code

Blogs - Hall of Fame

Free PDF Downloads

Quick Links

About DSPRelated.com

The Related Sites