Power Computation of a Digital Stream
clear all;
%example vector having ac & dc power
x = complex(randn(1,2048)+.2,randn(1,2048)-.1);
%total power, three equivalent methods
pwr1_total = mean(abs(x).^2); %mean of squared values
pwr2_total = mean(real(x).^2 + imag(x).^2);
pwr3_total = mean(x .* conj(x));
%total expressed as rms
rms_total = sqrt(pwr1_total);
%dc power
pwr1_dc = mean(real(x))^2 + mean(imag(x))^2; %square of mean of values
%ac power
pwr1_ac = pwr1_total - pwr1_dc; %mean of squares - square of mean
%relation of ac power to statistical variance
pwr2_ac = var(x); %approximately
%ac expressed as rms
rms1_ac = sqrt(pwr1_ac);
%ac relation to standard of deviation, std = sqrt(var)
rms2_ac = std(x); %approximately
%dc relation to variance
pwr2_dc = pwr1_total - var(x); %approximately
fprintf('----------------------------------------------------\r');
fprintf('power(total), : %1.5f, %1.5f, %1.5f\r',pwr1_total,pwr2_total,pwr3_total);
fprintf('rms(total) : %1.5f\r',rms_total);
fprintf('power(ac), variance : %1.5f, %1.5f\r',pwr1_ac,pwr2_ac);
fprintf('rms(ac), std : %1.5f, %1.5f\r',rms1_ac,rms2_ac);
fprintf('power(dc), (total-var) : %1.5f, %1.5f\r',pwr1_dc,pwr2_dc);
TI C281x Timer Configuration
**********************************************************************/
#include "DSP281x_Device.h"
/**********************************************************************
* Function: InitTimer()
*
* Description: Initialize the CPU Timer 0 at a rate of 2 Hz
**********************************************************************/
void InitTimer(void)
{
//CpuTimer0Regs.TIM.all = 0x0; //Timer0 Counter Register
CpuTimer0Regs.PRD.all = 0x047868C0; //Period Register
CpuTimer0Regs.TCR.all = 0x4830; //Control Register
/*
bit 15 0: CPU-Timer Interupt Flag
bit 14 1: Timer Interrupt Enable
bit 13-12 0's: Reserved
bit 11-10 10: CPU-Timer Emulation Modes
bit 9-6 0's: Reserved
bit 5 0: CPU-Timer Reload bit
bit 4 0: CPU-Timer Status Bit
bit 3-0 0's: Reserved
*/
CpuTimer0Regs.TPR.all = 0x0000; //Prescale Register
CpuTimer0Regs.TPRH.all = 0x0000; //Prescale Register High
PieCtrlRegs.PIEIER1.bit.INTx7 = 1; //Interruption enabled
IER |= 0x0001; //PIE group enabled
StartCpuTimer0();//Reset of the Timer
}
getsinc.m - TMX Transmultiplexer (MATLAB version)
function [sf] = getsinc(n_etapas, L)
% Experimentally obtained synchronization factors
if(n_etapas == 1)
sf = 2;
elseif (n_etapas == 2)
sf = L;
elseif L == 4
if(n_etapas == 3) % 8 branches
sf = 22;
elseif(n_etapas == 4) % 16 branches
sf = 36;
elseif(n_etapas == 5) % 32 branches
sf = 82;
end
elseif L == 6
if n_etapas == 3 % 8 branches
sf = 20;
elseif n_etapas == 4 % 16 branches
sf = 30;
elseif n_etapas == 5 % 32 branches
sf = 80;
end
elseif L == 8
if(n_etapas == 3) % 8 branches
sf = 26;
elseif(n_etapas == 4) % 16 branches
sf = 32;
end
elseif L == 10
if(n_etapas == 1)
sf = 2;
elseif(n_etapas == 2)
sf = 10;
end
end
PCM Encoding-Converting Quantized decimal sample values in to binary
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)
2D Image FIR Subsample
A = imread('in.jpg');
M = [ 1 4 6 4 1;
4 16 24 16 4;
6 24 36 24 6;
4 16 24 16 4;
1 4 6 4 1 ] / 256;
Y = filter2(M,A);
B = uint8(conv2(double(Y),M));
imwrite(B,'out.jpg','JPG');
Discrete Cosine Tranform Matrix generation
function[DCT] = dct_mtx(n)
[cc,rr] = meshgrid(0:n-1);
//disp(cc)
//disp(rr)
DCT = sqrt(2 / n) * cos(%pi * (2*cc + 1) .* rr / (2 * n));
DCT(1,:) = DCT(1,:) / sqrt(2);
endfunction
FIR BPF - Window Based
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//Band PASS FILTER - Window Based
clear all;
clc;
close;
M = 7 //Filter length = 7
Wc = [%pi/5,3*%pi/4]; //Digital Cutoff frequency
Wc2 = Wc(2)
Wc1 = Wc(1)
Tuo = (M-1)/2 //Center Value
hd = zeros(1,M);
W = zeros(1,M);
for n = 1:M
if (n == Tuo+1)
hd(n) = (Wc2-Wc1)/%pi;
else
n
hd(n) = (sin(Wc2*((n-1)-Tuo)) -sin(Wc1*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
if(abs(hd(n))<(0.00001))
hd(n)=0;
end
end
hd;
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp(h,'Filter Coefficients are')
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR BPF using Rectangular window M=7')
xgrid(1)
Average power of 1D signal
%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
Convolutional Encoding
//Caption:Convolutional Code Generation
//Time Domain Approach
close;
clc;
g1 = input('Enter the input Top Adder Sequence:=')
g2 = input('Enter the input Bottom Adder Sequence:=')
m = input('Enter the message sequence:=')
x1 = round(convol(g1,m));
x2 = round(convol(g2,m));
x1 = modulo(x1,2);
x2 = modulo(x2,2);
N = length(x1);
for i =1:length(x1)
x(i,:) =[x1(N-i+1),x2(N-i+1)];
end
x = string(x)
//Result
//Enter the input Top Adder Sequence:=[1,1,1]
//Enter the input Bottom Adder Sequence:=[1,0,1]
//Enter the message sequence:=[1,1,0,0,1]
//x =
//!1 1 !
//! !
//!1 0 !
//! !
//!1 1 !
//! !
//!1 1 !
//! !
//!0 1 !
//! !
//!0 1 !
//! !
//!1 1 !
Interpolation Demo
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
FFT at arbitrary frequencies (non-uniform FFT)
Following you find the code of a demo script (Demo.m) and that of a function (nufft.m). Just place them in separate m files (Demo.m and nufft.m) and execute the first.
------------------------------------------------------
Code of file Demo.m
------------------------------------------------------
echo on
%Author: J. Selva.
%Date: April 2011.
%Define a random vector with 1024 elements.
M = 1024;
a = rand(M,1)+1i*rand(M,1);
%Call the nufft function. The output is a struct with some variable for later use.
st = nufft(a)
%The field st.aF is the fft of size st.K.
%Define a vector of random frequencies.
f = rand([2*M,1])-0.5;
%Compute the DFT at frequencies in f and the derivatives of first and second order using
%loops and direct evaluation. This takes a while
%
tic; outDL = nufft(st,f,2,'directLoop'); tDL = toc;
%The function returns the DFT values in the first column, the first derivatives in the
%second column, and so on.
outDL(1:10,:)
%Do the same computation using interpolation. This is fast.
%
tic; outIL = nufft(st,f,2); tIL = toc;
%Compare the execution times
%
[tDL,tIL]
%The executions are faster if the code is vectorized. This is an example
%
tic; outDV = nufft(st,f,2,'directVec'); tDV = toc;
tic; outIV = nufft(st,f,2,'baryVec'); tIV = toc;
[tDV, tIV]
%For larger data sizes the difference is more significant. To check this, let us repeat
%the same experiment with M=1024*5.
%
M = 1024*5;
a = rand(M,1)+1i*rand(M,1);
st = nufft(a);
f = rand([2*M,1])-0.5;
tic; outDV = nufft(st,f,2,'directVec'); tDV = toc; %This line takes a while.
tic; outIV = nufft(st,f,2); tIV = toc;
%These are the timings:
[tDV,tIV]
%Finally, let us check the accuracy of the interpolated values
max(abs(outDV-outIV))./max(abs(outDV))
echo off
----------------------------------------------------------
----------------------------------------------------------
Code of file nufft.m.
----------------------------------------------------------
----------------------------------------------------------
function varargout = nufft(varargin)
%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, see
%
%J. Selva, 'Design of barycentric interpolator for uniform and nonuniform sampling
% grids', IEEE Trans. on Signal Processing, vol. 58, n. 3, pp. 1618-1627, March 2010.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
% st = nufft(a); %a is the vector of samples.
%
% The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
% 1/st.K.
%
% The DFT is defined using
%
% A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
% where m_1=-ceil(M/2) and M the number of samples.
%
% -Then call nufft with sintax
%
% out = nufft(st,f,nd,method)
%
% where
%
% f: Vector of frequencies where the DFT is evaluated. Its elements must follow
% abs(f(k))<=0.5
%
% nd: Derivative order. nufft computes derivatives up to order nd. At the output
% out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
% column of out contains the zero-order derivatives, i.e, the values of the DFT at
% frequencies in vector f.
% method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
% +'baryLoop': The output is interpolated using the barycentric method and a loop
% implementation.
%
% +'baryVec': The output is interpolated using the barycentric method and a
% vectorized implementation.
%
% +'directLoop': The output is computed using the DFT sum directly and a loop
% implementation.
%
% +'directVec': The output is computed using the DFT sum directly and a vectorized
% implementation.
if ~isstruct(varargin{1})
st = [];
st.a = varargin{1};
st.a = st.a(:);
st.M = length(st.a);
st.m1 = -ceil(st.M/2);
st.K = 2^ceil(log2(st.M)+1);
st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
%this would reduce st.P. The number of interpolation samples is 2*st.P+1.
st.T = 1/st.K;
st.B = -2*st.m1;
st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
st.vt = MidElementToEnd((-st.P:st.P)*st.T);
st.alpha = MidElementToEnd(baryWeights(st.T,st.B,st.P));
varargout = {st};
return
end
[st,f] = deal(varargin{1:2});
nd = 0;
if nargin > 2
nd = varargin{3};
end
method = 'baryVec';
if nargin > 3
method = varargin{4};
end
Nf = length(f);
out = zeros(Nf,nd+1);
switch method
case 'baryLoop' %Interpolated method using loops
for k = 1:length(f)
x = f(k);
n = floor(x/st.T+0.5);
u = x - n * st.T;
da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');
out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);
end
case 'baryVec' %Vectorized interpolated method
f = f(:);
Nf = length(f);
n = floor(f/st.T+0.5);
u = f - n * st.T;
pr = [-st.P:-1 , 1:st.P , 0];
ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
if length(f) == 1
ms = ms.';
end
out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
case 'directLoop' % Direct method using loops
for k = 1:length(f)
out(k,1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
end
for kd = 1:nd
out(k,kd+1) = 0;
for r = st.m1:st.m1+st.M-1
out(k,kd+1) = out(k,kd+1) + ...
((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
end
end
end
case 'directVec' %Vectorized direct method
for k = 1:length(f)
out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
for kd = 1:nd
out(k,kd+1) = ...
((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
end
end
end
varargout = {out};
function v = MidElementToEnd(v)
ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];
function v = APPulse(t,B,TSL)
v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));
function w = baryWeights(T,B,P)
vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;
N = length(vt);
LD = ones(1,N);
for k = 1:N
LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end
w = gam ./ LD;
w = w / max(abs(w));
function out = DerBarycentricInterp3(alpha,s,t,tau,nd)
vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,nd+1);
for k = 1:Nt-1
vD = vD*(tau-t(k))+alpha(k)*LF(k);
LF(k+1) = LF(k)*(tau-t(k));
end
vD = vD * (tau-t(Nt)) + alpha(Nt) * LF(Nt);
z = s;
for kd = 0:nd
z1 = z-z(end); cN = 0;
for k = 1:Nt-1
cN = cN * (tau-t(k))+z1(k)*alpha(k)*LF(k);
end
cN = cN/vD;
ztau = z(end)+(tau-t(end))*cN;
out(kd+1) = ztau;
if kd < nd
z = [ (z(1:end-1)-ztau)./(t(1:end-1)-tau) , cN ] * (kd+1);
end
end
function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)
NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory
%usage.
NBlock = 1+floor(NtauL/LBtau);
Nt = length(t);
out = zeros(NtauL,nd+1);
for r = 0:NBlock-1
ind1 = 1+r*LBtau;
ind2 = min([(r+1)*LBtau,NtauL]);
Ntau = ind2-ind1+1;
z = zL(ind1:ind2,:);
tau = tauL(ind1:ind2);
vD = zeros(Ntau,1);
LF = [1,zeros(1,Nt-1)];
LF = LF(ones(Ntau,1),:);
for k = 1:Nt-1
vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
LF(:,k+1) = LF(:,k) .* (tau-t(k));
end
vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);
for kd = 0:nd
pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);
for k = 1:Nt-1
cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
end
cN = cN ./ vD;
ztau = z(:,end)+(tau-t(end)) .* cN;
out(ind1:ind2,kd+1) = ztau;
if kd < nd
pr1 = ztau;
pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
pr2 = t(1:end-1);
pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
z = [pr1 ./ pr2, cN] * (kd+1);
end
end
end
Step Response Pole Effect Demo 2
for omega = 0.5:0.5:20
sys = tf([omega^2],[1,omega,omega^2]);
p = [1 omega omega^2];
r = roots(p);
hold on
subplot(1,2,1),plot(real(r), imag(r), 'o'), title('Pole Locations')
hold on
subplot(1,2,2),step(sys)
pause(0.5)
end
Step Response Pole Effect Demo
for zeta = 0.1:0.1:2
sys = tf([9],[1,6*zeta,9]);
p = [1 6*zeta 9];
r = roots(p);
hold on
subplot(1,2,1),plot(real(r), imag(r),'o'), title('Pole Locations'), xlabel('Real'), ylabel('Imaginary')
hold on
subplot(1,2,2),step(sys)
pause(.5)
end
WAV 2D Color Spectrogram
%FFT window size
window = 267;
[Y,FS,NBITS]=wavread('input.wav');
zz = zeros(0,window);
Y = Y.*1000;
for m = window+1:window:length(Y)
y = Y(m-window:m);
yy = fft(y);
zz = [zz, yy];
end
colormap(jet(280));
image(abs(zz));
xlabel('Time')
ylabel('Frequency')
DFT/IDFT Demo
clf;
L = 100;
for (L_factor = 1:1:6)
%L_factor = 3;
N_factor = 1;
L = round(L*L_factor);
n = [0:1:L-1];
x = sin(pi/1000*(n.*n));
subplot(3,1,1);
stem(n,x);
title ('x[n]');
subplot(3,1,2);
y = fft(x,L*N_factor);
plot(abs(y));
title ('DFT');
subplot(3,1,3);
xr = ifft(y);
stem(n,xr(1:L));
title ('IDFT');
pause;
end
Waterfall Spectrograph
function waterfallspect(s, fs, spectsize, spectnum)
%WATERFALLSPECT Display spectrum of signal s as a 3-D plot.
% s - input signal, use wavload to load a WAV file
% fs - Sampling frequency (Hz).
% spectsize - FFT window size
% spectnum - number of windows to analyze
frequencies = [0:fs/spectsize:fs/2];
offset = floor((length(s)-spectsize)/spectnum);
for i=0:(spectnum-1)
start = i*offset;
A = abs(fft(s((1+start):(start+spectsize))));
magnitude(:,(i+1)) = A(1:spectnum/2+1);
end
waterfall(frequencies, 0:(spectnum-1), magnitude');
xlabel('frequency');
ylabel('time');
zlabel('magnitude');
Image Quantizer Demo
clf;
for (truncate = 0:8)
x = imread('aqua.jpg');
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
x_t = (x/2^truncate) * 2^truncate;
%if above line doesn't work your MATLAB is old! use next line:
% x_t = uint8(double(uint8((double(x)/2^truncate))) * 2^truncate);
imshow(x_t);
subplot(2,2,3);
imshow(uint8(abs(double(x_t) - double(x))));
pause(0.5);
end
Signal-to-Distortion SDR Ratio
[f, fs, nbits] = wavread('input.wav');
[f2, fs2, nbits2] = wavread('input-truncated.wav');
top = 1/size(f,1)*(sum(f.^2))
bottom = 1/size(f2,1)*(sum(f2.^2))
Y = 10*log10(top/bottom)
RGB (JPEG) Image Separator
clf;
x = imread('input.jpg');
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);
subplot(2,2,1);
imshow(x);
title('Original');
subplot(2,2,2);
imshow(red);
title('Red');
subplot(2,2,3);
imshow(green);
title('Green');
subplot(2,2,4);
imshow(blue);
title('Blue');
Image Quantizer
%number of bits to truncate
tred = 4;
tgreen = 2;
tblue = 3;
clf;
x = imread('input.jpg');
red = x(:,:,1);
green = x(:,:,2);
blue = x(:,:,3);
red_t = (red/2^tred) * 2^tred;
green_t = (green/2^tgreen) * 2^tgreen;
blue_t = (blue/2^tblue) * 2^tblue;
x_t(:,:,1) = red_t;
x_t(:,:,2) = green_t;
x_t(:,:,3) = blue_t;
imshow(x);
pause;
imshow(x_t)
size = (24 - tred - tgreen - tblue) / 24