% *****************************************************************
% 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
#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;
}
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
% *****************************************************************
% 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'})
// 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))
function [y,X,H] = conv2d2(x,h)
[x1,x2] = size(x);
[h1,h2] = size(h);
X = zeros(x1+h1-1,x2+h2-1);
H = zeros(x1+h1-1,x2+h2-1);
Y = zeros(x1+h1-1,x2+h2-1);
for i = 1:x1
for j = 1:x2
X(i,j)=x(i,j);
end
end
for i =1:h1
for j = 1:h2
H(i,j)=h(i,j);
end
end
disp(X,'x=')
disp(H,'h =')
Y = fft2d(X).*fft2d(H);
disp(Y)
y = ifft2d(Y);
endfunction
/////////////////////////////////////////////
function [a2] = fft2d(a)
//a = any real or complex 2D matrix
//a2 = 2D-DFT of 2D matrix 'a'
m=size(a,1)
n=size(a,2)
// fourier transform along the rows
for i=1:n
a1(:,i)=exp(-2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a(:,i)
end
// fourier transform along the columns
for j=1:m
a2temp=exp(-2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a2(j,:)=a2temp.'
end
for i = 1:m
for j = 1:n
if((abs(real(a2(i,j)))<0.0001)&(abs(imag(a2(i,j)))<0.0001))
a2(i,j)=0;
elseif(abs(real(a2(i,j)))<0.0001)
a2(i,j)= 0+%i*imag(a2(i,j));
elseif(abs(imag(a2(i,j)))<0.0001)
a2(i,j)= real(a2(i,j))+0;
end
end
end
/////////////////////////////////////////////////
function [a] =ifft2d(a2)
//a2 = 2D-DFT of any real or complex 2D matrix
//a = 2D-IDFT of a2
m=size(a2,1)
n=size(a2,2)
//Inverse Fourier transform along the rows
for i=1:n
a1(:,i)=exp(2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a2(:,i)
end
//Inverse fourier transform along the columns
for j=1:m
atemp=exp(2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a(j,:)=atemp.'
end
a = a/(m*n)
a = real(a)
endfunction
function[SB_PSK] = PowerSpectra_MPSK()
rb = input('Enter the bit rate=');
Eb = input('Enter the energy of the bit=');
f = 0:1/100:rb;
Tb = 1/rb; //Bit duration
M = [2,4,8];
for j = 1:length(M)
for i= 1:length(f)
SB_PSK(j,i)=2*Eb*(sinc_new(f(i)*Tb*log2(M(j)))^2)*log2(M(j));
end
end
a=gca();
plot2d(f*Tb,SB_PSK(1,:)/(2*Eb))
plot2d(f*Tb,SB_PSK(2,:)/(2*Eb),2)
plot2d(f*Tb,SB_PSK(3,:)/(2*Eb),5)
xlabel('Normalized Frequency ---->')
ylabel('Normalized Power Spectral Density--->')
title('Power Spectra of M-ary signals for M =2,4,8')
legend(['M=2','M=4','M=8'])
xgrid(1)
endfunction
//Result
//Enter the bit rate in bits per second:2
//Enter the Energy of bit:1
//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)))
//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')
//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')
//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');
//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');
//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 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 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
