Details

Published: Tuesday, 17 September 2019 19:42

Hits: 3372

matlab: matlab is proprietary product of Mathworks. 1st version was released in 1984. It's full form is matrix laboratory, as it primarily dealt with matrix computations. MATLAB allows matrix manipulations, plotting of functions and data in 2D/3D, implementation of algorithms, interfacing with programs written in other languages, including C, Python, etc. Although MATLAB is intended primarily for numerical computing, an additional package Simulink, adds graphical multi-domain simulation. We learn matlab since we'll be using it a lot in signal processing, image/video processing, etc.

One of the best places to learn basic Matlab is from here: https://www.tutorialspoint.com/matlab/index.htm

Other useful resource is from matlab user uide or usinh help section of matlab product (click on help->Product_help on top bar after invoking matlab).

Invoking matlab: matlab needs to be installed on your m/c, before you can run it. Usually corporations/universities have licenses of Matlab, so you can use them there. Buying Matlab is expensive, isntead use open source equivalent.I'm presenting Matlab tutorial, since matlab syntax and feel is available in GNU Octave, which is open source equiv of Matlab. So, learning of Matlab can be used in Octave.

matlab -R2016a => starts matlab with version = R2016a. Every year matlab releases 2 versions a and b. Latest release is R2019a.
*.m file are matlab cmd files that can be run from cmdline in matlab. These are scripts that contain all matlab cmds.
ex: matlab_cmd1.m => just type matlab_cmd1 (without .m extension) on prompt of matlab
ex: matlab -nodisplay < cmd1.m => to run cmd1.m without opening gui

After starting matlab, we can start typing cmds on matlab command window. It has >> prompt for entering any cmds.

There is also a "workspace" window on top right, and "Command History" window on bottom right. Workspace window shows all var, while history window shows all cmds entered so far on prompt.

Matlab syntax:

spaces: Matlab is insensitive to spaces. It can parse code correctly even w/o spaces. spaces are not required b/w operators, but added for clarity.

comments: To comment out single line, use %. To comment out multiple lines use:
%{
 .... code ....
%}

format: format cmd changes precision of numbers. i.e format short;

Data types: Matlab is loosely typed, meaning we don't need to define dta type, matlab figures it out automatically.

1. Boolean: true(1) or false (0)

2. Numeric: Integer (signed/unsigned), floating point (single precision and double precision)

3. Text: char/string. string composed of array of char.

variables: A valid variable name starts with a letter, followed by letters, digits, or underscores. MATLAB is case sensitive, so A and a are not the same variable. In Matlab, all vars are arrays, and numeric var are of type double (i.e double precision).

>> a=1; => assigns value of 1 to var "a". If ";" not used at end of cmd, then results are displayed in cmd window (i.e a=1 will be displayed in cmd window). "a" is 1x1 array (array with 1 row and 1 column).  

>>5 => any cmd that returns an o/p val that is not assigned to a var, is assigned to reserved var "ans" (i.e answer). So, in this case, ans is assigned 5.

>> whos a => "whos" on any var shows it's attributes. here it shows "a" as 1x1 array of type "double". It shows 8 bytes are used to store this array "a" (since double precision requires 8 bytes). Just "whos" w/o any var will show info for all var.

>> B = [12, 62, 93, -8, 22; 16, 2, 87, 43, 91; -4, 17, -72, 95, 6] => creates 3x5 array (3 rows x 5 columns). Each row has 5 values. ; indicates end of row. After last row, no ; needed. [ ] indicates it's an array.

To get element from 2nd row, 3rd col of an array, we do B(2,3). NOTE: arrays are defined with square brackets, but accessed via round brackets. Curly brackets {} are 3rd kind of brackets used for creating cell array (explained later).

>> A=magic(5); => creates a 5x5 array with random entries in it.

>> C = [1,2,3] => an array with 1 row and 3 columns. C=[1;2;3] creates an array with 3 rows and 1 column.

>> Mystring = 'hello world'; => strings are char array of 1 row x n cols. single quotes used to store strings. Each char stored in 2 bytes, so this char array of 1x11 stored in 22 bytes.

>> A=ones(5)*4.7, B=magic(5); C=A/B; => multiple cmds can be entered on same line separated by , or ;. Commands that end with a comma display their results, while commands that end with a semicolon do not.

>> s= 1 -1/2 + ...

0.2 => gives result as 0.7. ellipses (...) allow continuation of stmt to next line

cell array: A cell array is a data type with indexed data containers called cells, where each cell can contain any type of data. Cell arrays commonly contain either lists of text, combinations of text and numbers, or numeric arrays of different sizes. Refer to sets of cells by enclosing indices in smooth parentheses (). Access the contents of cells by indexing with curly braces {}

ex: mycell = {'Mary', [3], 'summer', true}. To extract contents of mycell, we use {}. eg. mycell{1} gives Mary, mycell{2} or mycell{1,2} gives 3 and so on. When we use round brackets, we get cell array which is subset of cell array. i.e mycell(1) gives 'Mary' which is a 1x1 array, mycell(2) or mycell(1,2) gives [17] which is again 1x1 array, etc. NOTE: with curly braces, we got extracted values w/o any '..', or [..].

cell(3) => creates a 3x3 cell array of empty matrices.

cell arrays are very useful to store diff kind of data and then to access them in a loop to process them.

operators: 3 types of main operators:

1. Arithmetic operators: to perform numeric computations. +, -, .*(numeric multiplication), ./(numeric right division), .\(numeric left division), .^(numeric power), *(matrix multiplication), /(matrix right division), \(matrix left division), ^(matrix power), etc. Both operands must have same dimension. If one operand is scalar and the other is not, then MATLAB applies the scalar to every element of the other operand—this property is known as scalar expansion. right division (/ or ./) refers to regular division i.e if A=2, B=3, then A/B or A./B = 0.6667. left division (\ or .\) refers to inverted division - A\B or A.\B = B/A = 3/2 =1.5.

2. Relational operators: to compare operands. <, >, ==(equal to), ~=(not equal to), etc. Both operands must have same dimension. operation is applied to each element, and result is for each element, as o(false) or 1 (true). If one operand is scalar and the other is not, then MATLAB applies the scalar against every element of the other operand.

3. logical operators: these can be applied bit-wise, element wise or short circuit. &, |, ~, xor (i.e xor(A,B))

4. Assignment operator: = is an assignment operator as seen above.

5. colon operator (:) => most useful operator of matlab.

>>1:10 => creates a row vector, containing integers from 1 to 10, ans = [1 2 .... 10]

>>x=0:0.2:10 => creates an array with 1 row and 50 columns, with value of x going from 0 to 10 in crements of 0.2, i.e x=[0 0.2 0.4 0.8 .... 50]

: operator within round brackets used to extract specific rows/column entries from arrays. ex:

>> a(:,2) => extracts 2nd col of an array. This col retains it's structure, i.e final result is array with n rows and 1 column

>> a(3,:) => extracts 3rd row of an array. This row retains it's structure, i.e final result is array with 1 row and n columns

>> a(:,2:3) => extracts second and third column of array. final result is array with n rows and 2 columns

>> a(2:3,2:3) => extracts second and third rows and second and third columns. final result is array with 2 rows and 2 columns

Conditional Stmt: All conditional stmt end with keyword "end". 3 kinds:

1. if stmt => optional "elseif" "else" can also be used

ex: if a<30 disp 'small'; else disp 'large'; end => We should write it in multiple lines. If we write all of it in 1 line, we'll need ; to separate out multiple lines of "if"

2. case stmt: switch ..case ... end.

3. for loop => to loop specific num of times

ex: for n=2:6 x(n)=x(n-1); end => repeats with n=2,3,4,5,6

4. while loop => loops as long as condition is true. break (to exit) and continue (to skip to next iteration) stmt can be used from inside while loop.

ex: while fact<100 n=n+1;...; end

functions: can be called by typing name of func with appr args inside round brackets.

>> maxA = max(A,B); here max func is called with args A and B, both of which are arrays. o/p stored in var MaxA. If func returns more than 1 o/p, then we can store them by assigning them to separate var

>> helpFile = which('help');
>> [helpPath,name,ext] = fileparts(helpFile); => here fileparts func returns 3 o/p. We store them in 3 different var. To save memory space, we may want to ignore storing some o/p. We can do that by assigning ~ to that o/p. i.e [~,~,ext] = fileparts(helpFile);

Date functions => "date", "now", "calender", etc func display current date, etc.

import/export functions => to import/export data from external files.

1. importdata func: ex: filename = 'smile.jpg'; A = importdata(filename); image(A); => imports image file and then displays image in window

2. save func: ex: save my_data.out num_array -ascii => my_data.out is the delimited ASCII o/p data file created, num_array is i/p numeric array and −ascii is the specifier.

3. saveas: saveas(gcf,"plot_tmp"); saves graph generated as "plot_tmp".

misc functions:

1. ones, zeros: returns 1 or 0.

ex: y=ones; => returns scalar 1.

ex: y=ones(5); => returns 5x5 matrix with all entries as 1. We could also write it as y=ones(5,5);

ex: y=zeros(2,3) => returns 2x3 array of 0

print/display functions =>

1A. disp(x) => display value of any var (array, string, etc).
ex: A[10 20]; S='Hello'; name='alice', age=12; S=['my name', name,'age',num2str(age)]; => everything in S needs to be string
    disp(A); disp(S); disp('Hello'); => When arg of any function is a text string, we enclose it in single quotes ' ... '. disp 'hello' => is equally valid as we can use cmd syntax instead of func syntax.

1B. sprintf => used to create text. can be stored in a var.
ex: name='alice', age=12; S=sprintf('my name %s age %d', name,age);
    disp(S);

1C. fprintf => used to directly display the text without creating a variable. Always put \n at end. fscanf/fprintf works like C scanf/printf cmds.
ex: fprintf('my name %s age %d \n', name,age); => directly displays this without needing disp

Ploting: Ploting 2D/3D graphs is the most useful feature in Matlab, and used a lot.

1. plot: plot(x,y) plots y against x

ex: x = [-100:5:100]; y = x.^2; plot(x, y) => this plots 2D graph with y=x^2 on y axis and x from -100 to +100 (in increments of 5) on x axis. If increment reduced to 1, then graph becomes smoother as there are more points.

We can add title, labels along the x-axis and y-axis, grid lines and also to adjust the axes to spruce up the graph.

plot(x, y), xlabel('x'), ylabel('x^2'), title('sqaure Graph'), grid on, axis equal, axis ([-20, 20, 0,400]); => this plots x,y graph with label on x axis (as x), y axis (as y^2) and then a title on top of graph. It also adds grid lines on graph. "axis equal" causes equal spacing on both x and y axis (i.e linear plot). axis ([...]) sets min/max values for x,y axis. Here, xmin=-20, xmax=+20, ymin=0, ymax=400.

To plot multiple graphs on same plot, do plot(x1,y1, x2,y2, ...). ex: y = sin(x); g = cos(x); plot(x, y, x, g);

To have separate colors for each graph, add attr in plot(x,y, attr). For color attr, use r=red, g=green, b=blue, w=white, k=black, etc. ex: plot(x, y, 'r', x, g, 'g'); For diif line style, do: plot(x, y1, 'r', x, y2, 'b--o'); => creates lines as --o--o--o--

We can also have legends for graph (using legend cmd). We can print text on graph using text cmd

ex: text(10,-50, 'my custome text'); => prints the given text at x=10,y=-50 on graph. NOTE: x,y are coords in whatever unit is there on axis. If x is in freq(log scale) and y is in db, then x=10Hz,y=-50db is where the text printing starts.

ex: plot(y); => plots data in y verses index of each value. If y is a vector, then the x-axis scale ranges from 1 to length(y). ex:

subplot: With same figure, we can have various plots called subplots. We specify positions of various subplots using subplot cmd. First, specify subplot position for 1st graph, using subplot, then use plot(x,y) to plot 1st graph, then specify position for 2nd graph using subplot, then use plot(x,y) to plot 2nd graph, and so on ...

2. Bar: bar(x,y) plots bar chart for y against x.

3. Contour: A contour line of a function of two variables is a curve along which the function has a constant value. To draw contour map for given func g=f(x,y), first we need to create a meshgrid which defines a set of (x,y) points over domain of func. meshgrid defines x,y range over which we plot function g. Then, we use contour func: contour(x,y,g);

4. surf: Used to create 3D surface plot. Instead of contour, we draw 3D plot of func g=f(x,y). Here again we create meshgrid, then use surf func.

ex: [x,y] = meshgrid(-2:.2:2, -3:0.1:3); g = x .* exp(-x.^2 - y.^2); surf(x, y, g); => Here, func g is plotted on z axis. range of x is from -2 to +2, while y is -3 to +3.

5. ezplot: plot a given func, over range of values specified

ex: ezplot('x^2',[-1,2]) => plots square of x, over range of x=-1 to x=2.

ex: ezplot('x^2 - y^4') => plots x^2-y^4=0 over range x,y=-2Π to +2Π. This is the default range, when range is not specified.

Solving Equations: matlab can solve polynomial eqns using solve func:

ex: solve('(x-3)^2*(x-7)=0') => This solves the cubic eqn, and produces result as ans=[3; 3; 7] = 3x1 array. "=0" iin eqn s optional, as it's understood.

ex: eq = 'y^4 - 7*y^3 + 3*y^2 - 5*y + 9 = 0'; s = solve(eq,'y'); Now s becomes a 4x1 array. Each root can be accessed via s(0), s(1), s(2) and s(3), or as an array s. By default, eqn is solved for x, but to solve for other var, specify that var "y" in solve func.

We can also solve system of eqn.

ex: s = solve('5*x + 9*y = 5','3*x - 6*y = 4'); s.x; s.y; => displays soln for x,y

Calculus: Calc limit, solve differential/integral.

1. limit func: Calculates limit of func. The limit function falls in the realm of symbolic computing. Symbolic computing means using symbols as "x", etc in the o/p instead of getting numerical results. You need to use the syms function to tell MATLAB which symbolic variables you are using. You will get an error, if the function expects symbolic values, and you omit "syms" func for that var.

ex: syms x; limit((x^3 + 5)/(x^4 + 7),2); => 1st arg of limit func is the function, and 2nd arg is what the variable tends to i.e limit x-> a. Here we calc limit for this func, as x approaches 2. If 2nd arg not given, it's assumed to be 0.

2. differential: compute symbolic derivative using func diff().

ex: syms t; f = 3*t^2 + 2*t^(-2); diff(f,2); => this gives derivative as "6 - 12/t^2" (in symbols instead of numerical value). 2nd arg "2" in diff implies 2nd derivative (F''). If 2nd arg not provided, it's assumed to be 1 (1st derivative, F').

We can solve diff eqn symbolically using dsolve. ex: dsolve('D2y - y = 0','y(0) = -1','Dy(0) = 2')

3. integral: inverse of differentiation. "int" cmd calculates indefinite integral of any expression

ex: syms x; int(2*x); => this returns integral as x^2, though in reality integral is x^2+c, where c is arbitrary constant

ex: int(2*x,4,9) => This calculates definite integral with low limit=4, and hi limit=9. Since this yields numerical result, and not in symbols, we don't define x as symbol (i.e syms x; not used here)

Polynomial: To evaluate polynomial at any value, we use polyval func

ex: p = [1 7 0 -5 9]; polyval(p,4); Here p is defined as an array. When used inside polyval func, it represents polynomial P(x)=x^4+7*x^3-5*x+9. It's value is evaluated at x=4, which yields 693.

polyvalm evaluates matrix polynomial.

roots func calculates roots of poly. ex: p = [1 7 0 -5 9]; r = roots(p) => calc 4 roots of the poly eqn, and stores in var "r"

poly func => inverse of roots func, and returns polynomial coefficients

polyfit func => finds coeff of polynomial that fits a set of data in least square sense. We provide numerical values, x and y=f(x), and try to get poly func with max degree "n" that fits y over given values of x.

ex: x = [1 2 3 4 5 6]; y = [5.5 43.1 128 290.7 498.4 978.67]; p = polyfit(x,y,4); => here we ask matlab for poly with max degree "4". It returns 5 coeff for x^4, x^3, ...

Transform: Transform functions are ones that maps one domain to other domain. Depending on func used to do transform, it can be quite complex. matlab computes symbolic results.

1. Laplace transform: laplace transform converts f(t) in "t" domain into L{f(t)} in "s" domain. func laplace used for this.

ex: syms s t a w; laplace(a); laplace(t^2); laplace(sin(w*t)) => here we define s, t, a, w as symbols. "s" is for laplace transform o/p var, while others are i/p var going into laplace func. gives results as 1/s^2, 2/s^3 and w/(s^2+w^2)

2. inverse Laplace transform: ilaplace func used. Converts from s domain to t domain.

ex: syms s t ; ilaplace(s/(s^2+4)) => gives cos(2*t) as o/p

3. z transform:

4. Inverse z transform:

5. Fourier transform: converts from time domain to freq domain. freq domain expressed in radians(ω) or hertz(2Πf).

F=fourier(f,u,v); => calculates FT of f(u) as F(v), i.e F(v) = ∫ f(u) e^-jvu du. By default, if we omit u,v (i.e F=fourier(f)), then i/p func f is taken as f(x), and o/p func F is given as F(w), i.e F(w) = ∫ f(x) e^-jwx dx. Here, w is freq. F=fourier(f,w) imples f(x) and F(w). We have to declare u, v as syms (symbols), else fourier function will give an error as it expects symbols as it's inputs.

Generally, we use function syms fourier(f,t,w) which yields F(w)=  ∫ f(t) e^-jwt dt

 ex: syms x; f = exp(-2*x^2); FT = fourier(f); ezplot(FT); => shows FT = (2^(1/2)*pi^(1/2)*exp(-w^2/8))/2, and displays it on graph, with w=-2Π to +2Π.

6. Inverse Fourier transform: converts from freq domain to time domain. ifourier func used

ex: f = ifourier(-2*exp(-abs(w))); => gives result as f = -2/(pi*(x^2 + 1))

7. Fourier series: There is no function for calculating FS of periodic func. Instead, we directly compute terms C_k using integral in matlab and plot them. However, better way is to plot FT. Since FT is a superset of FS, we can look at FT values at distinct freq and get FS from that. Most of the times, FT suffices.

Generating basic signals: We can generate signals by using built in func. Basic code is:

x=[0:0.005:0.2]; % choose x axis points here
y=func(x); % generate y axis points by applying func to x
plot(x,y); % plot it

1. impulse function: y=dirac(x). Since dirac is a limiting func, whose value is 0 at x≠0, but infinity at x=0 (for continuous function), matlab draws dirac func with missing value at x=0, and 0 everywhere else.

2. sine wave: y=sin(x).

ex: create sine wave of 10Hz freq, T=0.1sec
t=[0:0.005:0.2]; => time is from 0s to 0.2s in steps of 5ms
Amp=1; freq=10; => Amplitude=1, freq=10Hz,
y=Amp*sin(2*pi*freq*t); => plots 2 full waves in 0.2sec, since 10Hz freq => T=0.1sec =10 sine wave/sec
plot(t,y); => plots t on x axis and sine o/p on y axis. plots over range of t from 0 to 0.2sec. plot is little choppy, but can be made smooth by choosing finer steps of 1ms, i.e t=[0:0.001:0.2]

If we choose freq=1000, then sine func argument=2*pi*1000*0.005*n=10*pi*n. So y=0 for all values of n. Then plot(t,y) should be 0 everywhere, since we don't have enough samples to sample any value other than 0. However, what we observe in reality is a choppy plot with y=very small value (10^-13 or so) and it bonces around 0. This is because sin(2*pi)=0 in theory, but computers calculate pi to finite digits only. That results in sin(2*pi) close to 0, but not exactly 0. Whenever you get weird plots like this, always check to make sure you have enough sample points.

ex: create sum of 2 sine waves of 10Hz freq, 

t=[0:0.005:0.2]; => time is from 0s to 0.2s in steps of 5ms
A=1;f1=10;f2=20 => freq=10Hz,20Hz
y=A*sin(2*pi*f1*t) + A*sin(2*pi*f2*t) ; => adds 2 sine waves
plot(t,y) => plots sum of 2 waves, amplitude of final "y" goes from -2 to +2

ex: show gibbs phenomemon in constructing square wave from sine waves (taken from mathworks website)

t = 0:.02:6.28; % t is 1Xn array
y = zeros(10,length(t)); %y is 10Xn array
x = zeros(size(t)); % x is 1Xn array
for k = 1:2:19
   x = x + sin(k*t)/k;
   y((k+1)/2,:) = x;
end
plot(y(1:2:9,:)')

------------------

1D. print text on graph
text_to_print = ['SNDR (in db) =',num2str(SNDR1),char(10),'SFDR (in db) =',num2str(SFDR1)]; => char(10) is needed to introduce newline. \n doesn't work.
text(10,-50,text_to_print); => prints the text above on graph.
saveas(gcf,"plot_tmp");

4. matlab native fft function are fft() and inverse is ifft(). To compute fft:
Y = fft(X) => computes the discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm. FFT is just bunch of points, it doesn't have any frequency info. We get this frequency info from sampling rate.
Y = fft(X,n) returns the n-point DFT. If no value is specified, Y is the same size as X. If X is a vector and the length of X is less than n, then X is padded with trailing zeros to length n. This is equiv to taking DFT of signal multiplied with a rectangular window. NOTE: for n-point dft, Y is also n point, i.e Y[0] to Y[n-1]. It's periodic with period n, so, only n points reported for Y.

4. PSD of signal
pxx = pwelch(x) => returns PSD of signal x using Welch's overlapped segment averaging estimator. If x is real-valued, pxx is a one-sided PSD estimate. If x is complex-valued, pxx is a two-sided PSD estimate.
[pxx,f] = pwelch(x,window,noverlap,f,fs) returns the two-sided Welch PSD estimates at the frequencies specified in the vector, f.

5. plotspec( ) => This function is defined as below (not std function of matlab). It computes the PSD of a signal, and plots Power in db on Y axis and freq on X axis (log scale).
 
plotspec.m script => uses pwelch signal and then draws plot
------
function [P,F]=PlotSpec(Signal,Fs,Navg,EnFig,color)
 FftLength=round(length(Signal)/Navg);
 Win=window(@blackmanharris,FftLength);

 if FftLength>length(Signal),
  warning('FftLength is greater than the signal length.')
 end

 [Ps,F] = pwelch( sin([0:FftLength-1]'*2*pi/16),Win, floor(3*FftLength/4), FftLength,Fs);
 norm=max(Ps(:,1));
 [P,F] = pwelch( Signal,Win, 0, FftLength,Fs); => note: P is not in db here.
 P=P(:,1)/norm;

if (EnFig==1)
    semilogx(F,10*log10(abs(P)),color,'LineWidth',2);
elseif (EnFig==2)
    plot(F,10*log10(abs(P)),color,'LineWidth',2);
end
    grid on;
    ylabel('Mag (dB)','FontSize',13);
    xlabel('Frequency (Hz)','FontSize',13);

end
-----
ex: [P,F]= PlotSpec(signal, Fs,1,1,'k'); => In John's script
ex: plotspec(signal,Fs,'on',[1 100]); => sampling freq=Fs, freq range on x axis = 1Hz to 100Hz. we can also specify limits separately in other cmd.
ex: [magnitudes, frequencies] = plotspec(signal,Fs,'off'); => instead of plotting, we can also store values in variables


----------------------
Fs1 = 8388608;
digin1 = dlmread('data1.csv'); => reads csv file, put 1 data in each row
figure('Color','w','Name','SDMout-Dec (SV)','WindowStyle','docked');
%PlotSpec(digin1,Fs2,1,1,'k'); => comment starts with %
PlotSpec(digin1,Fs1,1,1,'k');
ylim([-140 0]); => y db limit
xlim([0 Fs1/2]); => x freq limit
------------------------
ex: create sine wave
t=[0:0.005:0.2]; => time is from 0s to 0.2s in steps of 5ms
A=1;f1=10000;f2=50000 => freq=10kHz,50Khz
y=A*sin(2*pi*f1*t) + A*sin(2*pi*f2*t) ; => adds 2 sine waves
plot(t,y) => plots t on x axis and y on y axis

ex: plot fft of sine wave:
f0 = 307;    //sinusoid freq
Fs = 44100; //Sampling frequency in Hz:
Ts = 1/Fs; //Sampling time:
M = 1500; //1500 samples of sinusoid => think of this as the length of a rectangular window that multiplies our infinitely long sinusoid sequence.
n = 0:(M-1); //array of num from 0 to 1499
y = cos(2*pi*f0*Ts*n); //y[n]=cos[2*pi*f*t*n] => draws approx 10 full cosine waveforms = f0/Fs*1500=10. do "plot(n,y);" to see it
Y = fft(y); //computes M point fft since there are M points in y
f = (0:(M-1))*Fs/M; //compute freq points. from 0 to Fs in 1500 steps
plot(f,abs(Y),'*') //plot f on x-axis and mag of Y on y axis, with * as discrete points
xlim([0 Fs/2])
title('Magnitude of DFT')
xlabel('Frequency (Hz)')

-----------------------
ex: save plots in loop
sdm_files = {
         'adc_data_out_ram_raw_k4_b.txt',
         'adc_data_out_ram_raw_k5_b.txt',
         'adc_data_out_ram_raw_k6_b.txt'};

freq=8388608;
%freq=(8384, 789, 123); => makes it an array with integer values accessed thru freq{k}

for k = 1:3
   sdm_file_name = strcat(path, sdm_files{k});
   display(sprintf('Processing %s', sdm_file_name));
   sdm_output = dlmread(sdm_file_name);
   title = sprintf('SDM Ouput k=%d', k+3);
   figure('Color','w','Name', title,'WindowStyle','docked');
   PlotSpec(sdm_output, SDM_freq,1,1,'k');
   ylim([-160 0]);
   xlim([0 SDM_freq/2]);
   out_file_name = strcat(k, '_filt_fbw.jpg');
   saveas(gcf,out_file_name);
end

---------------------------------
end

Details

Published: Tuesday, 17 September 2019 18:20

Hits: 3559

Signal processing refers to processing any i/p signal to generate a transformed o/p signal. The black box that does this transformation may be analog or digital circuit. The i/p, o/p signals may be continuous or discrete.

In Signal theory, we refer to signals as varying in time, (i.e signal F(t)). However signals can be function of any var "x" (i.e signal F(x)). It so happens that for any signal of interest in real scenario, x happens to be time "t". Although "t" is real, it can be complex too. Thus signal F(x) can be real or complex, with x itself being real or complex. As you can guess, in real world life scenario, signal F(x) is real.

Joseph Fourier proved that any periodic continuous function can be represented exactly via infinite sum of discrete sine and cosine waveforms. This series came to be known as Fourier series (FS). For non-periodic continuous signals, instead of discrete freq components of sine and cosine, it was a continuous freq spectrum, and came to be known as Fourier Transform (FT). This was very important discovery of 18th/19th century by Joseph Fourier, who introduced it for purpose of solving heat equations in a metal plate.

These are some good links for learning Fourier:

http://www.thefouriertransform.com/ => Superb explanation of fourier

Few imp formula before we get into FS:

1. freq in radian (ω) = 2Π*F = 2Π/T, where F = freq in hertz, T=time period in sec. Ex: linear freq of 1 Hz = 1 rotation/sec = 2*pi radians/sec = 6.28 radians/sec. This is very important formula to remember, as many textbooks compute Fourier in F, while many do it in ω. There is a constant factor 2Π that appears depending on whether you are dealing in Hz or radians.

2. Euler's formula: e^ix = cos(x) + i sin(x) => This is regarded as most remarkable formula in mathematics, as it establishes relationship b/w trignometric functions and complex exponential functions. "i" is defined as square root of -1 ( i = √-1 ) In maths, we use "i", but in physics, we use "j". Euler's formula can be proved easily, by writing taylor's expansion of exponential, and then separating it in 2 series, with 1 series corresponding to sin(x), while other series corresponding to cos(x). But it was Euler who saw it first.

3. Euler identity is Euler's formula with x=Π: e^iπ = cos(Π) + i sin(Π) = -1. For x=0 or 2Π, e⁰ = e^i2π = cos(2Π) + i sin(2Π) = 1. From this, it implies that e^jω₀T₀ = e^i2π = 1.

4. Complex number: Any complex number c=a+ib can be represented in magnitude and phase form. a+ib = mag.e^iθ where mag = √(a² + b²), and tan(θ) = b/a. So, c=|c|.e^i∠c. This is how we'll represent complex numbers in signal theory.

5. Impulse function or delta dirac function: This is a hypothetical function which is 1 at t=0, but is 0 everywhere else. This can be thought of as a gaussian function, with limit t->0. i.e Limit t->0 ∫ x(t) dt =1, where x(t) is gaussian function. This impulse function is used extensively in signal theory, since mathematically impulse function allows conversion from cont domain to discrete domain. NOTE: x(t) is a continuous func. The resulting integral is the area of that func, and thus is a cont func of "t". As t->0, in continuous form, the value at t=0 approaches ∞ for impulse function, as base "x" axis time scale reduces to 0, to keep area = 1. in discrete form, impulse function is defined as x[n]=1 which is 1 at n=0, and 0 everywhere else.

Impulse function δ(x) defined as limit x->0 ∫ δ(x) dx =1

Just as we define impulse func at x=0, we can define it at any x=x₀. limit x->0 ∫ δ(x-x₀) dx =1. Here δ(x) is 1 at x=x₀, but 0 everywhere else. So, the gaussian func is just shifted from being centered at x=0, to being centered at x=x₀.

Few Sine/Cosine formula:

0. cos(x) = (e^ix + e^-ix )/2 ; sin(x) = (e^ix - e^-ix )/2i => derived by using Euler's formula above

1. sin(x) = cos(x-Π/2)

2. ∫ sin(2Πnt)*sin(2Πmt)dt = 0 for n≠m, and ∫ sin(2Πnt)*sin(2Πmt)dt = 1/2 for n=m. Here n, m are both integers, and limit of t is from 0 to 1

3. ∫ cos(2Πnt)*cos(2Πmt)dt = 0 for n≠m, and ∫ cos(2Πnt)*cos(2Πmt)dt = 1/2 for n=m. Here n, m are both integers, and limit of t is from 0 to 1

4. ∫ cos(2Πnt)*sin(2Πmt)dt = 0 for any n, m. Here n, m are both integers, and limit of t is from 0 to 1. These 3 eqn are used in proof of FS later.

Continuous time (CT) signals: here F(x) is a continuous function over all values of x, where -∞ < x < ∞. Fourier series or Fourier transform of CT signals is always aperiodic in freq domain.

1. Continuous Time Fourier Series (CTFS): For continuous and periodic signals in time domain, we can compute CTFS. This CTFS is discrete and aperiodic in freq domain.

Any cont periodic signal is represented as infinite summation of cos and sine functions. This is called FS representation of function. CTFS is the fundamental eqn, and all other series/transform for cont/discrete waveforms are derived from CTFS. CTFS is what was derived by Fourier, and is surprisingly very easy to prove.

x(t) = a₀ + ∑ a_k * cos(ω₀kt) + ∑ b_k * sin(ω₀kt) , where k = 1 to ∞. Here feq ω₀ = 2Π / T₀, where T₀ is the freq of periodic signal (1/T₀ = f₀). Here a₀ is the DC freq component, while ω₀ = fundamental or 1st harmonic, 2ω₀ = second harmonic and so on. Coefficients a₀, a_k, b_k are called Fourier coefficents, and are computed as follows. Integral is taken over 1 period of the periodic signal.

a₀ = 1/T₀ ∫ x(t) dt,

a_k = 1/T₀ ∫ x(t) cos(2Πkt/T₀) dt

b_k = 1/T₀ ∫ x(t) sin(2Πkt/T₀) dt

NOTE: This series contains only real terms, and is called real FS. k above is +ve, and goes from 1 to infinity. In reality, k = 0 to ∞. It just so happens, that for k=0, cos(2Πkt/T₀) =1, while sin(2Πkt/T₀) = 0. So, a₀ = 1/T₀ ∫ x(t) dt, while b₀ = 0.  No -ve k allowed. From now on, we'll refer to FS coeff as {a_k, b_k} with k from 0 to ∞. The motivation for a₀ is to provide dc offset, since all sine/cosine functions are always centered around X-axis with y=0 (i.e no dc offset, they start from 0 and go back to 0 over 1 cycle). In order to represent any offset from x-axis, we need a constant, which is what a₀ provides.

Also, x(t+T) = x(t) is true when calc thru FS, as ω₀T₀ = 2Π.

Proof: We can prove FS above very easily. Proof is presented very nicely here:  http://faculty.olin.edu/bstorey/Notes/Fourier.pdf.

Let's assume that FS is true, and we can indeed rep any function as sum of sine and cosine.

x(t) = a₀ + ∑ a_k * cos(ω₀kt) + ∑ b_k * sin(ω₀kt)

Multiply both sides of eqn by sin(2Πmt), and then integrating over limit 0 to 1, we get values of a₀, a_k and b_k. So that implies that eqn above is true since we can get coeff which satisfy this eqn.

We can make above series simpler. We put, a_k = d_k cos(θ_k), b_k = d_k sin(θ_k). We can do this since any 2 numbers on x and y axis can be represented on a circle with a certain radius and angle.

x(t) = ∑ d_k cos(θ_k) * cos(ω₀kt) + ∑ d_k sin(θ_k) * sin(ω₀kt)  = ∑ d_k cos(θ_k+ω₀kt)  where d_k = √(a_k² + b_k²), and tan(θ_k) = b_k /a_{k ,} with k = 0 to ∞.

Above eqn is easier to work with. It shows that any function can be rep as cosine series only with different magnitude and offset for each freq component. This seems logical, since sine can be converted to cosine with 90 degree phase offset. Instead of representing series in cosine series only, we could represent it in sine series only (as it's just phase shift of 90 degrees). Though above series, whether in sine or cosine, is ok to work in maths, we prefer to work in complex exponential in physics. This makes it easier to solve them. So, we try to get it in Euler's form.

We use the formula: cos(θ) = (e^jθ + e-^jθ)/2 to convert real numbers to complex numbers. We also note that subscript k can be both +ve/-ve value. Substituting this, we get the below form:

x(t) = ∑ C_k * e^jω₀kt ( -∞ < k < ∞ ),  Here C_k is complex number = d_k/2 * e^jθ_k = 1/2 * {d_k cos(θ_k) ₊ d_k j.sin(θ_k)}.

C_k = 1/T₀ ∫ x(t) e^-jω₀kt dt, Integral is taken over 1 period of the periodic signal.

Magnitude of this complex number C_k is d_k/2 = √(a_k² + b_k²)/2, while phase is θ_k where tan(θ_k) = b_k /a_{k ,}

NOTE: This series contains complex terms, and is called complex FS. Here k is both -ve and +ve, so we have freq spectrum on both sides of x axis (while FS terms were only on +ve x axis). This is why we have a 1/2 multiplier factor when going from real FS to complex FS, to account for it existing on both sides (we half the magnitude on both sides, when converting cosine to complex exponential). Also, FS coeff have -ve term in exponent, while inverse FS (i.e reconstruction of original signal) has +ve term in exponent. This is going to be true in all series and transforms as we'll see later (though it could have been other way around too, this is just a convention ?? FIXME?).

From now on, we will work with complex FS only. It's easy to represent, as we can draw magnitude and phase plot of this one complex number (instead of having 2 real numbers a_k and b_k). This plot runs over all values of k, so we see this plot on both sides of x axis. The magnitude plot of C_k is symmetrical around x=0 axis, while phase plot of C_k is -ve of each other on 2 sides of x-axis. Note that even though C_k is complex, the fourier series sum can still yield real function.

We can calculate FS coefficients for a variety of periodic signals, and observe many interesting things about sinusoidal wave, square wave, triangular wave, etc. We'll discuss those later.

2. Continuous Time Fourier Transform (CTFT): For continuous and aperiodic signals in time domain, we can compute CTFT. This CTFT is continuous and aperiodic in freq domain. Non periodic signal can be considered as a special case of a periodic signal, where Time period T₀ is expanded to infinity, and we consider only 1 cycle of periodic signal.

We can derive FS for aperiodic signals (attach proof, FIXME). On deriving, we come up with below eqn: NOTE: factor T₀ goes away and gets replaced by (2Π) when deriving this.

x(t) = 1/(2Π) ∫ X(ω) e^jωt dω, where X(ω) is the CTFT of x(t). Integral is taken over full -∞ < ω < ∞. Only difference compared to CTFS is that ω₀k is replaced by ω, and infinite summation is replaced by integral. This is also called inverse CTFT.

X(ω) = ∫ x(t) e^-jωt dt, Integral is taken over full time of the aperiodic signal -∞ < t < ∞. Note: this looks similar to C_k for periodic signals, except that C_k for aperiodic signals has become continuous now. Also, 1/T₀  is absent, as we folded that factor 1/(2Π) outside the integral in x(t) eqn above, so that it is no longer a part of FT. This is the forward CTFT. We can also write X(ω)  as X(jω) so as to be explicit that ω is complex.

So, to conclude, we see that CT periodic signals have discrete freq at +/-ω₀, +/-2ω₀, +/-3ω₀ and so on. But CT aperiodic signals have a continuum of freq from -∞ < ω < ∞. Other way to see is that freq is kω₀ for both periodic and aperiodic. But for periodic, k can take integer values only (i.e 0,1,2,3, ..), but for aperiodic signals, k can take all real values (i.e, 0, 0.001, 0.002, ....). CTFT is a superset of CTFS, as we can calculate CTFT of both periodic and aperiodic signals. There is nothing that prohibits taking CTFT of periodic signals. CTFT of aperiodic signals yields cont freq spectrum as we saw above. However, for periodic signals, we saw that we get discrete values for C_k. Since CTFT is derived from CTFS, CTFT should also yield same answer. However, we saw that CTFT is cont, and integral can never yield non cont discrete values. Dirac delta function comes to our rescue here, as they allow discrete values to be represented as the limiting case of gaussian func. So, CTFT of periodic signals will yield impulse or dirac delta at discrete freq, so that the func remains cont. So CTFT will serve us well for all kinds of waveforms (both periodic and aperiodic). Going forward, we'll learn CTFT only.

NOTE: we derived FT in terms of ω, which gives us a factor of (2Π). If we derive FT in terms of f, then we get rid of (2Π).

X(2Πf) = ∫ x(t) e^-j(2Πf)t dt => X(f) = ∫ x(t) e^-j(2Πf)t dt => remains same as X(ω), as integral is done over "t" and not "f"

x(t) = 1/(2Π) ∫ X(2Πf) e^j(2Πf)t d(2Πf) = 1/(2Π) * (2Π) ∫ X(2Πf) e^j(2Πf)t d(f) = ∫ X(2Πf) e^j(2Πf)t df => x(t) = ∫ X(f) e^j(2Πf)t df => gets rid of factor 2Π compared to X(ω), as integral is over "f".

3. Laplace transform: LT is defined for CT signals, in exactly same way as FT is defined for CT signals. In LT, complex var "s" has both real and imaginary parts (s = σ + jω). FT happens to a special case of LT, with σ=0, so s=jω. LT and FT are kind of same. When dealing with signals, we use FT, while when dealing with systems, we use LT. LT allows us to solve ordinary differential eqn encountered in RLC ckt, in an easier way by transforming it from time domain to laplace "s" domain (we can solve it equally easily by transforming it to Fourier "jω" domain). We eventually end up substituting s with jω, so both LT and FT end up same.

LT[x(t)] = X(s) = ∫ x(t) e^-st dt, Integral is taken over full time of the aperiodic signal -∞ < t < ∞. This is called two sided LT. Since most signals are 0 for t<0, we usually evaluate integral from 0 < t < ∞. This is called one sided LT. Note: this looks similar to FT eqn above, except that jω is replaced by s.

Discrete time (DT) signals: here F{x[n]} is a discrete function over all values of n, where x[n] is the input signal with n=0, +/-1, +/-2, ... ∞. Fourier series or Fourier transform of DT signals is always periodic in freq domain. This is the main diff b/w CT and DT. Otherwise DT are similar to CT.

1. Discrete Time Fourier Series (DTFS): For discrete and periodic signals in time domain, we can compute DTFS. This DTFS is discrete and periodic in freq domain. It's similar to it's CT counterpart CTFS.

We can derive DTFS from CTFS.

First let's derive x[n] from x[(t):

For CTFS, x(t) = ∑ C_k * e^jω₀kt ( -∞ < k < ∞ ). However, now x(t) is not continuous, but has discrete values. Assume N samples for x(t) equally spaced by T₀/N intervals, from t=0*T₀/N to t=(N-1)*T₀/N. So samples are x(0) to x((N-1)*T₀/N). 

x(t=n*T₀/N) = ∑ C_k * e^{j(2Π/T₀)k(n*T₀/N)} ( -∞ < k < ∞ ).

Represent discrete x(t) as x[n], where x[n] = x(n*T₀/N), and n = 0 to N-1.

x[n] = ∑ C_k * e^j(2Π/N)kn,  Here summation is over 1 period of C_k (Need to find out how summation went from ∞ to one period). Since C_k is periodic with period N (as we'll see below), we choose k as 0 < k < N-1. It doesn't matter what values of k we take, as long as it represents one full cycle over N values. So, C_k may have total N samples from C₀  to C_N-1 . i.e 0 < k < N-1, or C₁ to C_N , i.e 1 < k < N, or any other period. This looks similar to CTFS, except that it has period ω₀ replaced by 2Π/N. T₀ = 2Π / ω₀ = N.

Now, let's derive C_k for x[n] from C_k for x[(t):

For CTFS, C_k = 1/T₀ ∫ x(t) e^-jω₀kt dt, However, since x(t) is not continuous any more, but has discrete values, we can't calculate integral. However, we can make approximation of this integral, by taking only those "t" for which values of x(t) is known.

Then we use same value of x(t) for interval T₀/N. So, x(t) at any time t is x(n*T₀/N) where n = 0 to N-1.

C_k = 1/T₀ ∫ x(t) e^-jω₀kt dt ≈  1/T₀ ∑ x(n*T₀/N) e^{-jω₀kn*T₀/N}  * T₀/N, where t=n*T₀/N, ans summation is over n=0 to N-1. We approximate each integral piece as a rectangular box with height=x(n*T₀/N) and width=T₀/N.

=> C_k = 1/N ∑ x(n*T₀/N) e^{-j(2Π/T₀)kn*T₀/N} = 1/N  ∑ x[n] e^-j(2Π/N)kn, where x[n] = x(n*T₀/N), and n = 0 to N-1.

C_k = 1/N  ∑ x[n] e^-j(2Π/N)kn, summation is over 1 period of the periodic input signal with period=N samples, i.e 0 < n < N-1. C_k is complex number, and is periodic. with 0 < k < N-1. C_k is periodic can be deduced from the fact that for k=0 or k=N, C₀ = C_N as e^j(0) = e^j(2Π) . This looks similar to CTFS, except that it has summation instead of integral and C_k is periodic with period N.

Proof here: https://www.math.ubc.ca/~feldman/m267/dft.pdf

2. Discrete Time Fourier Transform (DTFT or DTF): For discrete and aperiodic signals in time domain, we can compute DTFT. This DTFT is continuous and periodic in freq domain. It's similar to it's CT counterpart CTFT.

We can derive DTFT from DTFS, the same way we derived CTFT from CTFS (by treating summation as integral when T -> ∞)

x[n] = 1/(2Π) ∫ X(ω) e^jωn dω, where X(ω) is the DTFT of x[n]. Integral is taken over 1 cycle of ω, i.e λ < ω < λ+2Π. Only difference compared to CTFT is that t is replaced by n, and integral over all values of ω, is replaced by integral over 1 cycle of ω. This is also called inverse DTFT.

X(ω) = ∑ x[n] e^-jωn, summation is taken over all samples of the aperiodic signal -∞ < n < ∞. Note: this looks similar to CTFT, except that t is replaced by n, and integral is replaced by summation, Also, X(ω) is periodic now, since X(ω+2Π) = X(ω). This is also called forward DTFT.

Since, same X(ω) notation is used for both CTFT and DTFT, it can cause confusion, as to whose FT is it representing (discrete or continuous). So, we for DTFT, we write it as X(Ω) or X(e^jω), while for CTFT, we write it as X(ω) or X(jω).

CTFT: X(ω) or X(jω) = |X(jω)|e^i∠X(jw).

DTFT: X(Ω) or X(e^jω) = |X(e^jω)|e^i∠X(ejω).

3. z transform (ZT): ZT is defined for DT signals, in exactly same way as LT is defined for CT signals. It is same as DTFT, except that complex var "z" is of form Ae^jω. DTFT is special case of ZT, where complex var "Ω" is of form e^jω (i.e z with unit magnitude only).

ZT{x[n]} = X(z) = ∑ x[n] z^-n,  -∞ < n < ∞. This is same as DTFT with z=e^jω.

ZT is used to solve discrete diff eqn in time domain. Concept is same as that of LT for solving cont diff eqn in time domain.

----

Discrete Fourier Transform (DFT):

This is a very popular topic in any signals theory. When we talk about anything in DSP (digital signal processing), DFT is what we are talking about. As we saw above, we can take DTFT of discrete signals which are not periodic. Since in real life, most signals are not periodic and are sampled, DTFT is what we use mostly for discrete signals (instead of DTFS). Computers can work with discrete signals (samples of analog signals at discrete times), and process them to generate DTFT of those samples. However, DTFT that we get above is itself continuous, and that is a problem for machines that deal with discrete samples only. Computers can't generate analog DTFT waveform. They can generate samples of DTFT. We can generate samples of DTFT so close by that we can approximate analog DTFT signal. Those samples of DTFT are called DFT.

DTFT of x[n] = X(Ω) = ∑ x[n] e^-jωn,

DFT of x[n] = X_k or X[k] = samples of X(Ω) at discrete freq =  X(Ω) | _ω=2Πk/N where 0 <= k <= N-1. Here we are taking N samples, X₀ to X_N-1

But discrete freq components implies periodic x[n], since only periodic signals can have discrete freq components. As we saw above, DTFS has discrete freq components, while DTFT has continuous freq spectrum. Basically as signal becomes aperiodic (or T -> ∞), more and more freq components start showing up, until it becomes continuous freq spectrum.

So, if we try to reconstruct original signal x[n] by taking inverse DFT of X_{k ,} we'll not get original x[n], but periodic x[n]. As we have more and more samples of X_k, we'll get more and more original values of x[n] with larger period (i.e T starts going to ∞). So, DFT of x[n] is basically same as taking DTFS of periodic x[n].

DTFS of x[n] = C_k = 1/N  ∑ x[n] e^-j(2Π/N)kn

DFT of x[n] = X_k or X[k] = 1/N  ∑ x[n] e^-j(2Π/N)kn = 1/N  ∑ x[n] e^-jω_kn , where ω_k =2Πk/N where 0 <= k <= N-1. 

Thus DFT of x[n] can be calculated by either of 2 ways:

1. Compute DTFT of x[n]. Then take samples of DTFT at discrete freq  ω_k =2Πk/N where 0 <= k <= N-1. 

2. Compute DTFS of x[n]. This DTFS gives same result as 1 (taking DTFT samples) above.

Inverse DFT of X_k = x[n] = ∑ X_k e^jω_kn, where ω_k =2Πk/N where 0 <= k <= N-1. 

This is same as inverse DTFS = x[n] = ∑ C_k * e^j(2Π/N)kn,  where 0 <= k <= N-1. 

Thus from now on, we'll treat DFT the same as DTFS. So, when we say DFT, we mean DTFS, and not DTFT.

Details

Published: Monday, 16 September 2019 22:21

Hits: 2942

vlsi digital standards

There are many standards for diffeent files used in vlsi. We talk about some standards and their formats:

LEF:

DEF:

UVM:

...

Details

Published: Monday, 16 September 2019 22:08

Hits: 2389

Webserver: Any website that you access, goes to a computer on internet. On that computer is a program running (for linux it's apache) which allows that computer to return the response back to the client. This computer is called a webserver. Before we get into webserver program, let's get some basics.

IP addr:

Each device connected to internet is assigned a unique 32 bit IP addr. These are written as 4 octets(each octet is 8 bits) - e.g. 18.251.48.77. IP addresses can go from 0.0.0.0 to 255.255.255.255. However, IP adr are divided into 5 classes based on addr range:

class A: 1.0.0.0 - 127.255.255.255 => Top 8 bits are Net ID. Bottom 24 bits are Node ID.

class B: 128.0.0.0 - 191.255.255.255 => Top 16 bits are Net ID. Bottom 16 bits are Node ID.

class C: 192.0.0.0 - 223.255.255.255 => Top 24 bits are Net ID. Bottom 8 bits are Node ID.

class D: 224.0.0.0 - 239.255.255.255 => multicast. Top 4 MSB are "1110". Remaining 28 bits are Multicast ID.

class E: 240.0.0.0 - 255.255.255.255 => experimental, reserved for future use

Public vs Private: IP addr can be public or private.