function [p]= RaisedCosineSpectrum()
//Practical Solution for Intersymbol Interference
//Raised Cosine Spectrum
rb = input('Enter the bit rate:');
Tb =1/rb;
t =-3:1/100:3;
Bo = rb/2;
Alpha =0; //roll-off factor Intialized to zero
x =t/Tb;
for j =1:3
for i =1:length(t)
if((j==3)&((t(i)==0.5)|(t(i)==-0.5)))
p(j,i) = sinc_new(2*Bo*t(i));
else
num = sinc_new(2*Bo*t(i))*cos(2*%pi*Alpha*Bo*t(i));
den = 1-16*(Alpha^2)*(Bo^2)*(t(i)^2)+0.01;
p(j,i)= num/den;
end
end
Alpha = Alpha+0.5;
end
a =gca();
plot2d(t,p(1,:))
plot2d(t,p(2,:))
poly1= a.children(1).children(1);
poly1.foreground=2;
plot2d(t,p(3,:))
poly2= a.children(1).children(1);
po1y2.foreground=4;
poly2.line_style = 3;
xlabel('t/Tb------>');
ylabel('p(t)------->');
title('RAISED COSINE SPECTRUM - Practical Solution for ISI')
legend(['ROlloff Factor =0','ROlloff Factor =0.5','ROlloff Factor =1'])
xgrid(1)
endfunction
//Result
//Enter the bit rate:1
/* -------------- SinCos.h begins -------------- */
#ifndef __SINCOS_H__
#define __SINCOS_H__
#include "types.h"
/*---------------------------------------------------------------------------;
; This function computes the sine of an angle using the Maclaurin series ;
; expansion method. Input angle is signed fractional where -1 corresponds to ;
; -pi and +1 to +pi; output is signed fractional and the minimum negative ;
; value is limited to -1+2^(-15) (0x8001 hexadeximal) to avoid overflow in ;
; subsequent calculations. ;
; ;
; Input: Y0 = angle (signed fractional; -1 = -pi, +1 = +pi) ;
; ;
; Output: Y0 = sine (signed fractional) ;
; ;
; Registers modified: A, B, X0, Y ;
;---------------------------------------------------------------------------*/
asm Frac16 Sin(Frac16 angle);
/*---------------------------------------------------------------------------;
; This function computes the cosine of an angle using the Maclaurin series ;
; expansion method. Input angle is signed fractional where -1 corresponds to ;
; -pi and +1 to +pi; output is signed fractional and the minimum negative ;
; value is limited to -1+2^(-15) (0x8001 hexadeximal) to avoid overflow in ;
; subsequent calculations. ;
; ;
; Input: Y0 = angle (signed fractional; -1 = -pi, +1 = +pi) ;
; ;
; Output: Y0 = cosine (signed fractional) ;
; ;
; Registers modified: A, B, X0, Y ;
;---------------------------------------------------------------------------*/
asm Frac16 Cos(Frac16 angle);
#endif //ifndef __SINCOS_H__
/* -------------- SinCos.h ends -------------- */
/* -------------- SinCos.c begins -------------- */
/*---------------------------------------------------------------------------;
; This function computes the sine of an angle using the Maclaurin series ;
; expansion method. Input angle is signed fractional where -1 corresponds to ;
; -pi and +1 to +pi; output is signed fractional and the minimum negative ;
; value is limited to -1+2^(-15) (0x8001 hexadeximal) to avoid overflow in ;
; subsequent calculations. ;
; ;
; Input: Y0 = angle (signed fractional; -1 = -pi, +1 = +pi) ;
; ;
; Output: Y0 = sine (signed fractional) ;
; ;
; Registers modified: A, B, X0, Y ;
;---------------------------------------------------------------------------*/
asm Frac16 Sin(Frac16 angle){
CLR.W Y1 //Clear Y1
MOVE.W Y0,B //Compute absolute value of angle
ABS B
CMP.W #$4000,B //|angle| > pi/2?
BLT SinAngleOK //No, proceed
ADD.W #$8000,Y0 //Yes, add pi and set flag for
BFSET #$0001,Y1 //final negation of result
SinAngleOK: MPYR Y0,Y0,A //Compute angle squared
MOVE.W A,X0 //X0 = angle squared
MOVE.W #$0002,B //Compute Maclaurin expansion
MOVE.W #$FFE2,A
MACR X0,B1,A //A1 = first partial result
MOVE.W #$0150,B
MACR X0,A1,B //B1 = second partial result
MOVE.W #$F669,A
MACR X0,B1,A //A1 = third partial result
MOVE.W #$28CD,B
MACR X0,A1,B //B1 = fourth partial result
MOVE.W #$AD52,A
MACR X0,B1,A //A1 = fifth partial result
MOVE.W #$3244,B
MACR X0,A1,B //B1 = sixth partial result
MOVE.W #$0003,X0 //Required shift amount
MPYR B1,Y0,A //A = result / 8
CMP #$1000,A //If magnitude is maximum
BNE SinCheckSatL //saturate the result
MOVE.W #$7FFF,A
BRA SinValOK
SinCheckSatL: CMP.W #$F000,A
BNE SinNoSat
MOVE.W #$8001,A
BRA SinValOK
SinNoSat: ASLL.L X0,A //Shift to get correct result
//(Maclaurin coefficients are
//divided by 8 to avoid overflow)
SinValOK: BRCLR #$0001,Y1,Update_Sin //Negate result if |angle| was
NEG A //greater than pi/2
Update_Sin: MOVE.W A,Y0 //Save result
RTS //Return from subroutine
}
/*---------------------------------------------------------------------------;
; This function computes the cosine of an angle using the Maclaurin series ;
; expansion method. Input angle is signed fractional where -1 corresponds to ;
; -pi and +1 to +pi; output is signed fractional and the minimum negative ;
; value is limited to -1+2^(-15) (0x8001 hexadeximal) to avoid overflow in ;
; subsequent calculations. ;
; ;
; Input: Y0 = angle (signed fractional; -1 = -pi, +1 = +pi) ;
; ;
; Output: Y0 = cosine (signed fractional) ;
; ;
; Registers modified: A, B, X0, Y ;
;---------------------------------------------------------------------------*/
asm Frac16 Cos(Frac16 angle){
ADD.W #$4000,Y0 //Add pi/2 to angle
JSR Sin //Call sine function
RTS //Return from subroutine
}
/* -------------- SinCos.c ends -------------- */
/* -------------- Usage example begins ------------- */
/* application specific includes */
#include "SinCos.h"
/* global variables */
Frac16 Angle,cosphi,sinphi;
/* Function calls */
cosphi=Cos(Angle); //Compute sin and cos of angle
sinphi=Sin(Angle); //(for use in Park transforms)
/* -------------- Usage example ends ------------- */
% 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)';
#
# rfft.py
#
# This file contains a recursive version of the fast-fourier transform and
# support test functions. This module utilizes the numpy (numpy.scipy.org)
# library.
#
# References
# - http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
# - "A Simple and Efficient FFT Implementation in C++", by Vlodymyr Myrnyy
import numpy
from numpy.fft import fft
from numpy import sin, cos, pi, ones, zeros, arange, r_, sqrt, mean
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def rFFT(x):
"""
Recursive FFT implementation.
References
-- http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
-- "A Simple and Efficient FFT Implementation in C++"
by Vlodymyr Myrnyy
"""
n = len(x)
if (n == 1):
return x
w = getTwiddle(n)
m = n/2;
X = ones(m, float)*1j
Y = ones(m, float)*1j
for k in range(m):
X[k] = x[2*k]
Y[k] = x[2*k + 1]
X = rFFT(X)
Y = rFFT(Y)
F = ones(n, float)*1j
for k in range(n):
i = (k%m)
F[k] = X[i] + w[k] * Y[i]
return F
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def getTwiddle(NFFT=8):
"""Generate the twiddle factors"""
W = r_[[1.0 + 1.0j]*NFFT]
for k in range(NFFT):
W[k] = cos(2.0*pi*k/NFFT) - 1.0j*sin(2.0*pi*k/NFFT)
return W
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def DFT(x, N=8):
"""
Use the direct definition of DFT for verification
"""
y = [1.0 + 1.0j]*N
y = r_[y]
for n in range(N):
wsum = 0 + 0j;
for k in range(N):
wsum = wsum + (cos(2*pi*k*n/N) - (1.0j * sin(2*pi*k*n/N)))*x[k]
y[n] = wsum
return y
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def test_rfft(N = 64, # FFT order to test
nStart = 0.2, # Note aliased signal included
nStep = 2.1, # Samples per period step
pStep = pi/4, # Phase step size
limErr = 10e-12, # Error limit to check
maxErr = 0 # Max difference
):
"""
Use the built in numpy FFT functions and the direct
implemenation of the DFT to verify the recursive FFT.
This testbench verifies the different implementations are within
a certain limit. Because of the different implemenations the values
could be slightly off (computer representation calculation error).
"""
# Use test signal nStart:nStep:N samples per cycle
for s in arange(nStart, N+nStep, nStep):
for p in arange(0, pi+pStep, pStep):
n = arange(N, 0, -1)
x = cos(2*pi*n/s + p)
xDFT = DFT(x,N)
nFFT = fft(x,N)
xFFT = rFFT(x)
rmsErrD = sqrt(mean(abs(xDFT - xFFT))**2)
rmsErrN = sqrt(mean(abs(nFFT - xFFT))**2)
if rmsErrD > limErr or rmsErrN > limErr:
print s, p, "Error!", rmsErrD, rmsErrN
print xDFT
print nFFT
print xFFT
if rmsErrD > maxErr:
maxErr = rmsErrD
elif rmsErrN > maxErr:
maxErr = rmsErrN
print "N %d maxErr = %f " % (N,maxErr)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If the module is run test a bunch of different size FFTs
if __name__ == '__main__':
# The following is fairly exhaustive and will take some time
# to run.
tv = 2**arange(1,12)
for nfft in tv:
test_rfft(N=nfft)
/* -------------- MovAvg.h begins -------------- */
#ifndef __MOVAVG_H__
#define __MOVAVG_H__
#include "types.h"
typedef struct
{
unsigned int decimation;
unsigned int deccnt;
unsigned int BufHeadPtr;
unsigned int BufLength;
Frac16 multiplier;
Frac16 ykm1;
}MovAvgparams;
/*---------------------------------------------------------------------------;
; Initialize buffer length, pointer and decimation divider of moving average ;
; filter; initialize memory buffer of moving average filter to all zeroes. ;
; ;
; Input: Y0 = buffer length ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = buffer length (filter order) ;
; R3 = pointer to storage buffer for the past input samples ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void MovAvgInit(unsigned int length, MovAvgparams *params, Frac16 *data);
/*---------------------------------------------------------------------------;
; This function implements a moving average Nth order FIR filter. ;
; The result yk is computed as the sum of the current and the past N input ;
; samples multiplied by the multiplier value. Canonical moving average ;
; filtering is obtained if multiplier = 1/(N+1). ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample) ;
; (R2) = decimation (1 = no decimation, 2 = discard every other ;
; sample, and so on) ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = N (length of buffer and filter order) ;
; (R2+4) = fractional multiplier ;
; (R2+5) = yk-1 (past output) ;
; R3 = Pointer to storage buffer for the N past input samples ;
; ;
; Output: Y0 = yk ;
; ;
; Registers modified: A, X0, Y, R0, R3 ;
;---------------------------------------------------------------------------*/
asm Frac16 MovAvg(Frac16 xk, MovAvgparams *params, Frac16 *data);
#endif //ifndef __MOVAVG_H__
/* -------------- MovAvg.h ends -------------- */
/* -------------- MovAvg.c begins -------------- */
/*---------------------------------------------------------------------------;
; Initialize buffer length, pointer and decimation divider of moving average ;
; filter; initialize memory buffer of moving average filter to all zeroes. ;
; ;
; Input: Y0 = buffer length ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = buffer length (filter order) ;
; R3 = pointer to storage buffer for the past input samples ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void MovAvgInit(unsigned int length, MovAvgparams *params, Frac16 *data){
MOVE.W #0001,X:(R2+1) //Initialize decimation divider
CLR.W X:(R2+2) //Initialize buffer head pointer
MOVE.W Y0,X:(R2+3) //Initialize buffer length
REP Y0 //Initialize past data to zero
CLR.W X:(R3)+
RTS //Return from subroutine
}
/*---------------------------------------------------------------------------;
; This function implements a moving average Nth order FIR filter. ;
; The result yk is computed as the sum of the current and the past N input ;
; samples multiplied by the multiplier value. Canonical moving average ;
; filtering is obtained if multiplier = 1/(N+1). ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample) ;
; (R2) = decimation (1 = no decimation, 2 = discard every other ;
; sample, and so on) ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = N (length of buffer and filter order) ;
; (R2+4) = fractional multiplier ;
; (R2+5) = yk-1 (past output) ;
; R3 = Pointer to storage buffer for the N past input samples ;
; ;
; Output: Y0 = yk ;
; ;
; Registers modified: A, X0, Y, R2, R3 ;
;---------------------------------------------------------------------------*/
asm Frac16 MovAvg(Frac16 xk, MovAvgparams *params, Frac16 *data){
DEC.W X:(R2+1) //Decrement decimation counter
BEQ DoMovAvg //Decremented to zero?
MOVE.W X:(R2+5),Y0 //No, output old value and exit
BRA MovAvgDone
DoMovAvg: MOVE.W X:(R2)+,X0 //Yes, reload decimation counter
MOVE.W X0,X:(R2)+
ADDA #1,SP //Preserve N (needed if called
MOVE.W N,X:(SP) //from an ISR)
MOVE.W X:(R2)+,A //Load index of buffer head into
MOVE.W A1,N //address modifier register
INC.W A X:(R2)+,Y1 //Point to next buffer location
//and load buffer length
CMP.W Y1,A //If buffer head index = buffer
BNE MAIndexOK //length, then reset buffer head
CLR.W A1 //index
MAIndexOK: MOVE.W A1,X:(R2-2) //Save buffer head index
MOVE.W X:(R3+N),A1 //Load oldest buffer sample
MOVE.W Y0,X:(R3+N) //Save new input sample
MOVE.W X:(R2)+,Y0 //Load multiplier
MPY A1,Y0,A X:(R3)+,X0 //Accumulate oldest sample
REP Y1 //Accumulate the latest N samples
MAC X0,Y0,A X:(R3)+,X0
MOVE.W A,Y0 //Return accumulation result
MOVE.W Y0,X:(R2) //Save new output sample
MOVE.W X:(SP)-,X0 //Restore N
MOVE.W X0,N
MovAvgDone: RTS //Return from subroutine
}
/* -------------- MovAvg.c ends -------------- */
/* -------------- Usage example begins ------------- */
/* application specific constants */
#define MovAvgOrder 40 //Order of moving average filter
#define MovAvgDecimation 100 //Decimation factor for m.a.
/* application specific includes */
#include "MovAvg.h"
/* global variables */
MovAvgparams MAparams;
Frac16 MAdata[MovAvgOrder];
int FiltIn,MAOut;
/* initializations */
MAparams.decimation=MovAvgDecimation;
MAparams.multiplier=FRAC16((float)1/(MovAvgOrder+1));
MovAvgInit(MovAvgOrder,&MAparams,(Frac16*)&MAdata); //Initialize m.a. filter
/* Moving average filter function call */
MAOut=MovAvg(FiltIn,&MAparams,(Frac16*)&MAdata); //Call moving av. filter
/* -------------- Usage example ends ------------- */
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
//Generation of Differential Phase shift keying signal
clc;
bk = [1,0,1,1,0,1,1,1];//input digital sequence
for i = 1:length(bk)
if(bk(i)==1)
bk_not(i) =~1;
else
bk_not(i)= 1;
end
end
dk_1(1) = 1&bk(1); //initial value of differential encoded sequence
dk_1_not(1)=0&bk_not(1);
dk(1) = xor(dk_1(1),dk_1_not(1))//first bit of dpsk encoder
for i=2:length(bk)
dk_1(i) = dk(i-1);
dk_1_not(i) = ~dk(i-1);
dk(i) = xor((dk_1(i)&bk(i)),(dk_1_not(i)&bk_not(i)));
end
for i =1:length(dk)
if(dk(i)==1)
dk_radians(i)=0;
elseif(dk(i)==0)
dk_radians(i)=%pi;
end
end
disp(bk,'(bk)')
bk_not = bk_not';
disp(bk_not,'(bk_not)')
dk = dk';
disp(dk,'Differentially encoded sequence (dk)')
dk_radians = dk_radians';
disp(dk_radians,'Transmitted phase in radians')
function [c] = PCM_Encoding(x,L,en_code)
//Encoding: Converting Quantized decimal sample values in to binary
//x = input sequence
//L = number of qunatization levels
//en_code = normalized input sequence
n = log2(L);
c = zeros(length(x),n);
for i = 1:length(x)
for j = n:-1:0
if(fix(en_code(i)/(2^j))==1)
c(i,(n-j)) =1;
en_code(i) = en_code(i)-2^j;
end
end
end
disp(c)
//PN sequence generation
//Maximum-length sequence generator
//Program to generate Maximum Length Pseudo Noise Sequence
clc;
//Assign Initial value for PN generator
x0= 1;
x1= 0;
x2 =0;
x3 =0;
N = input('Enter the period of the signal')
for i =1:N
x3 =x2;
x2 =x1;
x1 = x0;
x0 =xor(x1,x3);
disp(i,'The PN sequence at step')
x = [x1 x2 x3];
disp(x,'x=')
end
m = [7,8,9,10,11,12,13,17,19];
N = 2^m-1;
disp('Table Range of PN Sequence lengths')
disp('_________________________________________________________')
disp('Length of shift register (m)')
disp(m)
disp('PN sequence Length (N)')
disp(N)
disp('_________________________________________________________')
// Hamming Weight and Hamming Distance
//H(7,4)
//Code Word Length = 7, Message Word length = 4, Parity bits =3
close;
clc;
//Getting Code Words
code1 = input('Enter the first code word');
code2 = input('Enter the second code word');
Hamming_Distance = 0;
for i = 1:length(code1)
Hamming_Distance =Hamming_Distance+xor(code1(i),code2(i));
end
disp(Hamming_Distance,'Hamming Distance')
//Result
//Enter the first code word [0,1,1,1,0,0,1]
//Enter the second code word[1,1,0,0,1,0,1]
//Hamming Distance 4.
