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
//Caption: Reading a Speech Signal &
//[1]. Write it in another file
//[2]. Changing the bit depth from 16 bits to 8 bits
clear;
clc;
[y,Fs,bits_y]=wavread("E:\4.wav");
wavwrite(y,Fs,8,"E:\4_8bit.wav");
[x,Fs,bits_x]=wavread("E:\4_8bit.wav");
Ny = length(y); //Number of samples in y (4.wav)
Nx = length(x); //Number of samples in x (4_8bit.wav)
Memory_y = Ny*bits_y; //memory requirement for 4.wav in bits
Memory_x = Nx*bits_x; //memory requirement for 4_8bit.wav in bits
disp(Memory_y,'memory requirement for 4.wav in bits =')
disp(Memory_x,'memory requirement for 4_8bit.wav in bits =')
//Result
//memory requirement for 4.wav in bits =
// 133760.
//
//memory requirement for 4_8bit.wav in bits =
// 66880.
% By David Valencia
% As posted in DSPRelated.com
% at http://www.dsprelated.com/showcode/45.php
% Used for the DWT chain processing program
function [hx] = formfiltersdwt(n_stages,branch,h0,h1)
if(branch==1)
hx=formfilters(n_stages,branch,h0,h1);
elseif(branch==2)
hx=formfilters(n_stages,2,h0,h1);
else
hx=formfilters(n_stages-branch+2,2,h0,h1);
end
% The code line that is changed in formfilters.m for DWT is:
%p = mod(2^(n_stages-1-i),p);
/* -------------- PI32.h begins -------------- */
#ifndef __PI32_H__
#define __PI32_H__
#include "types.h"
typedef struct
{
Frac16 c0;
Frac16 c1;
Frac16 IntSatH;
Frac16 IntSatL;
Frac16 kp;
UInt16 nkp;
Frac16 OutSatH;
Frac16 OutSatL;
}PI32params;
typedef struct
{
Frac32 Skm1;
Frac16 xkm1;
}PI32data;
/*---------------------------------------------------------------------------;
; Initialize integral memory and past input of PI regulator. ;
; ;
; Input: A10 = Integral memory initialization value ;
; Y0 = Past input initialization value ;
; (R2) = Integral memory storage destination ;
; (R2+2) = Past input storage destination ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void PI32Init(Frac32 S, Frac16 xp, PI32data *data);
/*---------------------------------------------------------------------------;
; This function implements a 32 bit precision PI regulator with separate ;
; saturations on the integral portion and on the whole output. ;
; The routine first computes the integral portion at current sampling ;
; instant k (Sk) as Sk = Sk-1 + c0*xk + c1*xk-1. Sk is subsequently ;
; saturated between IntSatH and IntSatL and saved for later use at the next ;
; sampling instant. ;
; Next, proportional portion Pk is computed as kp*2^nkp*xk and saturated ;
; between -1 and +1. The two portions are added together and the result ;
; yk=Sk+Pk is saturated between OutSatH and OutSatL. ;
; Coefficients c0 and c1 are equal and amount to ki*Ts/2, where Ts is the ;
; sampling period and ki is the integral multiplier. The proportional ;
; multiplier kp*2^nkp is split into a fractional part and an exponent in ;
; order to allow values greater than unity. ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample of error signal) ;
; (R2) = c0 ;
; (R2+1) = c1 ;
; (R2+2) = IntSatH ;
; (R2+3) = IntSatL ;
; (R2+4) = kp ;
; (R2+5) = nkp ;
; (R2+6) = OutSatH ;
; (R2+7) = OutSatL ;
; (R3) = Sk-1 (32 bits) ;
; (R3+2) = xk-1 ;
; ;
; Output: A = yk (32 significant bits) ;
; ;
; Registers modified: A, B, X0, Y, R0 ;
;---------------------------------------------------------------------------*/
asm Frac32 PI32(Frac16 xk, PI32params *params, PI32data *data);
#endif //ifndef __PI32_H__
/* -------------- PI32.h ends -------------- */
/* -------------- PI32.c begins -------------- */
/*---------------------------------------------------------------------------;
; Initialize integral memory and past input of PI regulator. ;
; ;
; Input: A10 = Integral memory initialization value ;
; Y0 = Past input initialization value ;
; (R2) = Integral memory storage destination ;
; (R2+2) = Past input storage destination ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void PI32Init(Frac32 S, Frac16 xp, PI32data *data){
MOVE.L A10,X:(R2)+ //Initialize integral memory
MOVE.W Y0,X:(R2) //Initialize past input
RTS //Return from subroutine
}
/*---------------------------------------------------------------------------;
; This function implements a 32 bit precision PI regulator with separate ;
; saturations on the integral portion and on the whole output. ;
; The routine first computes the integral portion at current sampling ;
; instant k (Sk) as Sk = Sk-1 + c0*xk + c1*xk-1. Sk is subsequently ;
; saturated between IntSatH and IntSatL and saved for later use at the next ;
; sampling instant. ;
; Next, proportional portion Pk is computed as kp*2^nkp*xk and saturated ;
; between -1 and +1. The two portions are added together and the result ;
; yk=Sk+Pk is saturated between OutSatH and OutSatL. ;
; Coefficients c0 and c1 are equal and amount to ki*Ts/2, where Ts is the ;
; sampling period and ki is the integral multiplier. The proportional ;
; multiplier kp*2^nkp is split into a fractional part and an exponent in ;
; order to allow values greater than unity. ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample of error signal) ;
; (R2) = c0 ;
; (R2+1) = c1 ;
; (R2+2) = IntSatH ;
; (R2+3) = IntSatL ;
; (R2+4) = kp ;
; (R2+5) = nkp ;
; (R2+6) = OutSatH ;
; (R2+7) = OutSatL ;
; (R3) = Sk-1 (32 bits) ;
; (R3+2) = xk-1 ;
; ;
; Output: A = yk (32 significant bits) ;
; ;
; Registers modified: A, B, X0, Y, R0 ;
;---------------------------------------------------------------------------*/
asm Frac32 PI32(Frac16 xk, PI32params *params, PI32data *data){
MOVE.W R2,X0
MOVE.W X0,R0
MOVE.L X:(R3)+,A //Load Sk-1
MOVE.W X:(R0)+,X0 //Load c0
MOVE.W Y0,Y1 //Save xk
MAC Y0,X0,A X:(R0)+,Y0 X:(R3)+,X0 //c0*xk, c1, xk-1
MAC Y0,X0,A X:(R0)+,Y0 X:(R3)-,X0 //c1*xk-1, IntSatH, decr. R3
MOVE.W Y0,B //Load IntSatH into B
CMP A,B Y1,X:(R3)+ //Check upper satur. and save xk
BGT CheckSatLo //If Sk<IntSatH no need for sat.
MOVE.W Y0,A //else saturate and increment R0
ADDA #1,R0 //to skip IntSatL in the
BRA IntSatOK //parameter structure
CheckSatLo: MOVE.W X:(R0)+,B //Load IntSatL
CMP A,B //Check lower saturation
BLE IntSatOK //If Sk>IntSatL no need for sat.
MOVE.W B1,A //else saturate
IntSatOK: ADDA #-3,R3,R3 //Point to Sk storage location
MOVE.L A10,X:(R3) //Save Sk
MOVE.W X:(R0)+,Y0 //Load kp
MPY Y0,Y1,B X:(R0)+,Y0 //kp*xk, nkp
CLB B,X0 //Number of leading bits of
CMP Y0,X0 //result greater or equal than
NOP //desired shift count?
BLT PropSat //No, saturate
ASLL.L Y0,B //Yes, shift
BRA PropSatOK
PropSat: TST B //Saturate proportional portion
BLT PropSatNeg //according to sign
MOVE.W #$7FFF,B
BRA PropSatOK
PropSatNeg: MOVE.W #$8000,B
PropSatOK: ADD B,A X:(R0)+,B //Compute yk, load OutSatH
CMP A,B X:(R0)+,X0 //Check upper sat., load OutSatL
BGT CheckSatOutLo //If yk<OutSatH no need for sat.
MOVE.W B1,A //else saturate
BRA OutSatOK
CheckSatOutLo: CMP X0,A //Check lower saturation
BGT OutSatOK //If yk>OutSatL no need for sat.
MOVE.W X0,A //else saturate
OutSatOK: RTS //Return from subroutine
}
/* -------------- PI32.c ends -------------- */
/* -------------- Usage example begins ------------- */
/* application specific constants */
#define PIc0 100
#define PIc1 100
#define PIIntSatH 15000
#define PIIntSatL -15000
#define PIkp 0x4000
#define PInkp 1
#define PIOutSatH 30000
#define PIOutSatL -30000
/* application specific includes */
#include "PI32.h"
/* global variables */
PI32params PIparams;
PI32data PIdata;
int PIIn,PIOut;
/* initializations */
PIparams.c0=PIc0; //Initialize PI regulator
PIparams.c1=PIc1; //parameters and past data
PIparams.IntSatH=PIIntSatH;
PIparams.IntSatL=PIIntSatL;
PIparams.kp=PIkp;
PIparams.nkp=PInkp;
PIparams.OutSatH=PIOutSatH;
PIparams.OutSatL=PIOutSatL;
PI32Init(0,0, &PIdata);
/* PI32 regulator function call */
PIOut=(int)(PI32(PIIn,&PIparams,&PIdata)>>16); //Call PI regulator
/* -------------- Usage example ends ------------- */
function [m, b] = ls_linear(x_observed, y_observed)
%
% Perform Least Squares curve-fitting techniques to solve the coefficients
% of the linear equation: y = m*x + b. Input and Output equation
% observations must be fed into the algorithm and a best fit equation will
% be calculated.
%
% Usage: [M,B] = ls_linear(observed_input, observed_output);
%
% observed_input is the x-vector observed data
% observed_output is the y-vector observed data
%
% Author: sparafucile17 06/04/03
%
if(length(x_observed) ~= length(y_observed))
error('ERROR: Both X and Y vectors must be the same size!');
end
%Calculate vector length once
vector_length = length(x_observed);
%
% Theta = [ x1 1 ]
% [ x2 1 ]
% [ x3 1 ]
% [ x4 1 ]
% [ ... ... ]
%
Theta = ones(vector_length, 2);
Theta(1:end, 1) = x_observed;
%
% Theta_0 = [ y1 ]
% [ y2 ]
% [ y3 ]
% [ y4 ]
% [ ... ]
%
Theta_0 = y_observed;
%
% Theta_LS = [ M ]
% [ B ]
%
Theta_LS = ((Theta' * Theta)^-1 ) * Theta' * Theta_0;
%Return the answer
m = Theta_LS(1);
b = Theta_LS(2);
function save_fig(h, name, format)
% Matlab has a wierd way of saving figures and keeping the plot looking
% exactly the way it does on your screen. Thi function simplifies the
% entire process in a way that is barely documented.
%
% Usage: SAVE_FIG(h, name, format);
%
% H is the handle to the figure. This can be obtain in the
% following manner: H = figure(1);
% NAME is the filename without extension
% FORMAT is the graphic format. Options are:
%
% 'bmpmono' BMP monochrome BMP
% 'bmp16m' BMP 24-bit BMP
% 'bmp256' BMP 8-bit (256-color)
% 'bmp' BMP 24-bit
% 'meta' EMF
% 'eps' EPS black and white
% 'epsc' EPS color
% 'eps2' EPS Level 2 black and white
% 'epsc2' EPS Level 2 color
% 'fig' FIG Matlab figure file
% 'hdf' HDF 24-bit
% 'ill' ILL (Adobe Illustrator)
% 'jpeg' JPEG 24-bit
% 'pbm' PBM (plain format) 1-bit
% 'pbmraw' PBM (raw format) 1-bit
% 'pcxmono' PCX 1-bit
% 'pcx24b' PCX 24-bit color PCX file format
% 'pcx256' PCX 8-bit newer color (256-color)
% 'pcx16' PCX 4-bit older color (16-color)
% 'pdf' PDF Color PDF file format
% 'pgm' PGM Portable Graymap (plain format)
% 'pgmraw' PGM Portable Graymap (raw format)
% 'png' PNG 24-bit
% 'ppm' PPM Portable Pixmap (plain format)
% 'ppmraw' PPM Portable Pixmap (raw format)
% 'svg' SVG Scalable Vector Graphics
% 'tiff' TIFF 24-bit
%
% Author: sparafucile17
if(strcmp(format,'fig'))
saveas(h, name, 'fig');
else
options.Format = format;
hgexport(h, name, options);
end
function Demo
close all
T = 1; %Sampling period.
B = 0.7; %Two-sided bandwidth.
%Plot error spectrum for P = 10, 20, and 30.
PlotErrorSpectrum(T,B,10);
figure
PlotErrorSpectrum(T,B,20);
figure
PlotErrorSpectrum(T,B,30);
function PlotErrorSpectrum(T,B,P)
% Frequency grid specification in figure.
f0 = 0;
M = (2*P+1)*10;
Df = B/(2*M);
f = f0+Df*(0:M-1);
ESinc = zeros(size(f));
Eg = zeros(size(f));
%Look for the maximum error among different shifts u for each frequency in f.
for u = -0.5:1/((2*P+1)*5):0.5
HRef = exp(1i*2*pi*f*u); %Ideal spectrum.
HSinc = chirpzDFT(sinc((-P:P)+u/T),-P*T,T,f0,Df,M); %Sinc spectrum.
Hg = chirpzDFT(gPulseCoefs(u,T,B,P),-P*T,T,f0,Df,M); %g-pulse spectrum.
ESinc = max([ESinc;abs(HRef-HSinc)]); %Compute current max. error for sinc.
Eg = max([Eg;abs(HRef-Hg)]); %Compute current max. error for g.
end
plot(f,20*log10([ESinc;Eg].'))
xlabel('Freq. (f*T)')
ylabel('Maximum error (dB)')
legend('sinc pulse','g pulse','Location','NorthOutside')
pr = get(gca,'YTick');
text(0,2*pr(end)-pr(end-1),['B*T = ',num2str(B*T),', P = ',num2str(P)])
title('Error spectrum (maximum for all possible shifts u)')
grid on
function c = gPulseCoefs(u,T,B,P)
t = (-P:P)*T+u;
c = sinc(t/T).*real(sinc((1-B*T)*sqrt((t/T).^2-P^2)))/real(sinc(1i*(1-B*T)*P));
function X = chirpzDFT(x,to,Dt,fo,Df,M)
%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
% N
% X(k) = sum x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
% n=1
%
%Input parameters:
%
% x: Signal samples.
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
% M: Desired number of output samples.
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%
x = x(:).';
N = numel(x);
P = 2^ceil(log2(M+N-1));
%The next three lines do not depend on the vector x, and so they can be pre-computed if
%the time and frequency grids do not change, (i.e, for computing the transform of
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and
%three if the grids change.
a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse.
%Weigh x using a and perform FFT convolution with phi.
X = ifft( fft(x.*a,P) .* phi );
%Truncate the convolution tails and weigh using b.
X = X(N:M+N-1) .* b;
function xdft=fdft(x)
N=length(x);
xdft(1,N)=0; % defining a 1xN sequence with all the elements equals to zero.
%*************THIS IS THE DFT IMPLEMENTATION**************
%For each value of k, the Summatory function is computed over the interval 0<n<N-1.
for k=0:N-1,
for n=0:N-1,
xdft(k+1)=xdft(k+1) + x(n+1)*exp(-j*2*pi*k*n/N);
end
end
%NOTE: You may skip the use of for-loops by defining k and n as two N-points sequences. This is:
% k=0:N-1; n=0:N-1;
%The DFT is now computed in matrix form (as a matrix multiplication of the N-points Fourier Matrix and the discrete-time sequence x).
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]);