DSPRelated.com
   Code Snippets

WAV Butterworth Subsample

Ron ● April 1, 2011 ● Coded in Matlab 
[X, fs, nbits] = wavread('file.wav');
[N, Wn] = buttord(0.12, 0.125, 3, 60);
[B,A] = butter(N,Wn,'low');
Y = filter(B,A,X);
Z = downsample(Y,8);
wavwrite (Z, fs/8, nbits, 'file-sampled.wav');

Aliasing Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
% First plot:
subplot(2,2,1)
t = 0:0.005:10;
xa = 2*t.*exp(-t);
plot(t,xa);grid
xlabel('Time, sec');ylabel('Amplitude');
title('Continuous-time signal x_{a}(t)');

% Second plot:
subplot(2,2,2)
wa = 0:10/511:10;
ha = freqs(2,[1 2 1],wa);
plot(wa/(2*pi),abs(ha));grid;
xlabel('Frequency, Hz');ylabel('Magnitude');
title('|X_{a}(j\Omega)|');
axis([0 5/pi 0 2]);

% Third plot:
subplot(2,2,3)
T = .5;
n = 0:T:10;
xs = 2*n.*exp(-n);
k = 0:length(n)-1;
stem(k,xs);grid;
xlabel('Time index n');ylabel('Amplitude');
title('Discrete-time signal x[n]');

% Fourth plot:
subplot(2,2,4)
wd = 0:pi/255:4*pi;
hd = freqz(xs,1,wd);
plot(wd/(2*T*pi), T*abs(hd));grid;
xlabel('Frequency, Hz');ylabel('Magnitude');
title('|X_{\delta}(j\Omega)|');
axis([0 4/(2*T) 0 2])

Sampling Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
% First plot
t = 0:0.0005:1;
f = 5;
x_t = cos(2*pi*f*t);
subplot (4,1,1);
plot (t,x_t);
grid;
xlabel ('Time [sec]');
ylabel ('Amplitude');
title ('Continuous-time signal x(t)');
axis([0 1 -1.2 1.2])

% Second plot
n = 0:0.125:1;
x_n = cos(2*pi*f*n);
subplot(4,1,2);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.125 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs < bound] (T=0.125)');
axis ([0 1 -1.2 1.2]);

% Third plot
n = 0:0.1:1;
x_n = cos(2*pi*f*n);
subplot(4,1,3);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.1 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs = bound] (T=0.1)');
axis ([0 1 -1.2 1.2]);

% Fourth plot
n = 0:0.075:1;
x_n = cos(2*pi*f*n);
subplot(4,1,4);
y = zeros(1,length(t));
for i = 1:length(n)
    y = y + x_n(i)*sinc(t/0.075 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs > bound] (T=0.075)');
axis ([0 1 -1.2 1.2]);

Interpolation Demo

Ron ● April 1, 2011 ● Coded in Matlab 
clf;
for steps = 3:20
    x = 0:2*pi/steps:2*pi;
    y = sin(x);
    plot(x,y,'o-');
    grid;
    axis([0 2*pi -1.2 1.2]);
    pause(0.5)
end

2D Matrix Downsample

Ron ● April 1, 2011 ● Coded in Matlab 
function [ Y ] = downsample2d(M)
x = 2;
y = 2;

N = downsample(M,x);
N = N';
P = downsample(N,y);
P = P'
end

Decimate WAV file

Ron ● April 1, 2011 ● Coded in Matlab 
factor = 8;
[X, fs, nbits] = wavread('file.wav');
Y = decimate(X,factor);
wavwrite (Y, fs/factor, nbits, 'file-decimate.wav');

Frequency Sweep WAV Creator

Ron ● March 25, 2011 ● Coded in Matlab 
function [] = SineSweep(fstart, fstop, length, method, fs, bits, PHI)
%SINESWEEP Creates a custom sine wave sweep.
%   SineSweep(fstart, fstop, length, [method], [fs], [bits], [PHI])
%   fstart (Hz) - instantaneous frequency at time 0
%   fstop (Hz) - instantaneous frequency at time length
%   length (s) - the length of time to perform the sweep
%   method - 'quadratic', 'logarithmic', or 'linear' (default)
%   fs (Hz) - the sampling frequency (default = 48KHz)
%   bits - bit depth of each sample 8, 16 (default), 24, or 32
%   PHI (deg) - initial phase angle. cos = 0, sin = 270 (default)
%
%   Written by: Ron Elbaz
%   Date: February 15, 2011

%check and set missing parameters
if exist('method','var') ~= 1
    method = 'linear';
end
if exist('fs','var') ~= 1
    fs = 48000;
end
if exist('bits','var') ~= 1
    bits = 16;
end
if bits ~= 8 && bits ~= 16 && bits ~= 24 && bits ~= 32
    bits = 16;
end
if exist('PHI','var') ~= 1
    PHI = 270;
end

%make sure start and stop frequencies are valid
if fstop < fstart
    temp = fstart;
    fstart = fstop;
    fstop = temp;
end

%avoid aliasing
if fstart > fs/2
    fstart = fs/2;
end
if fstop > fs/2
    fstop = fs/2;
end

%create time vector
t = 0:(1/fs):length;

%scale to avoid clipping due to rounding error
Y = 0.9999*chirp(t,fstart,length,fstop, method, PHI);

%play the sound
% sound(Y,fs);

%store sweep to wav file
wavwrite(Y,fs,bits,'SineSweep.wav');

end

NLMS Adaptive filter

Senthilkumar ● March 24, 2011●1 comment ● Coded in Matlab 
%Normalized Least Mean Square Adaptive Filter
%for Echo Cancellation
function [e,h,t] = adapt_filt_nlms(x,B,h,delta,l,l1)
%x = the signal from the speaker 
%B = signal through echo path
%h = impulse response of adaptive filter
%l = length of the signal x
%l1 = length of the adaptive filter
%delta  = step size
%e = error
%h = adaptive filter impulse response
tic;
for k = 1:150
    for n= 1:l
        xcap = x((n+l1-1):-1:(n+l1-1)-l1+1);
        yout(n) = h*xcap';
        e(n) = B(n) - yout(n);
        xnorm = 0.000001+(xcap*xcap');
        h = h+((delta*e(n))*(xcap/xnorm));
    end
    eold = 0.0;
    for i = 1:l
        mesquare = eold+(e(i)^2);
        eold = mesquare;
    end
    if mesquare <=0.0001
        break;
    end
end
t = toc; %to get the time elapsed

Average power of 1D signal

Senthilkumar ● March 24, 2011 ● Coded in Matlab 
%Function to calculate the Average of a Voice signal
function [y] = pow_1(x)
N =  length(x);  % Length of voice signal 'x'
xold =0.0;       %initialize it to zero
for n = 1:N
    sumx = xold+(x(n)^2);
    xold = sumx;
end
y = sumx/N

Frequency Estimate using FFT

● February 2, 2011●1 comment ● Coded in Matlab 
%For a real-valued, single tone sine or cosine wave, estimate the frequency 
%using an fft. Resolution will be limited by the sampling rate of the tone
%and the number of points in the fft.

%Procedure:
%1. Generate the test tone with a given Fs and N
%2. Take the fft and keep the first half
%3. Detect the maximum peak and its index
%4. Equate this index to a frequency

%Error Estimate:
% (Fs/N)= Hz/bin

%Things to keep in mind:
%Adding noise to the input will skew the results
%windowing will help with the fft

%% Clear and close everything
clc;
close all;
clear all; %remove this line if you are trying to use breakpoints

%% Just copy into m file, save and run
Fs = 1e6; %1MHz
fi = 1.333e3; %1kHz
t = 0:1/Fs:.5;
y = sin(2*pi*fi*t);

%Plot the time and frequency domain. Be sure to zoom in to see the waveform
%and spectrum
figure;
subplot(2,1,1);
temp = (2/fi)*Fs;
plot(y);
xlabel('time (s)');
subplot(2,1,2);
sX=fft(y);
N=length(sX);
sXM = abs(sX(1:floor(end/2))).^2; %take the magnitude and only keep 0:Fs/2
plot(0:Fs/N:(Fs/2)-1,sXM);
xlabel('Frequency')
ylabel('Magnitude')

[vv, ii]=max(sXM); %find the index of the largest value in the spectrum

freqEst = (ii-1)*Fs/N; %units are Hz
resMin = Fs/N; %units are Hz

display(freqEst);%display to the command window the results
display(resMin);