Thursday 11 December 2014

How to fit (curve-fit) a histogram using Gaussian / Normal distribution, in Matlab?

1. Use 'histfit'


histfit(DATA,NBINS)


%By default it will fit the histogram with the normal distribution.


2. There are lots of other distributions supported by the 'histfit'. The details are given in this Matlab link. Quick look of supported distributions from this link:

 3. Example of 'histfit' with normal distribution for GPS residuals.








Thursday 4 December 2014

How to have all the functions in a single file, instead of having individual .m file for each function, in Matlab?

1. It is easy. It can be done by using functions-even the top level script should also be written as a function. An example is given below:

All these functions are present in Top.m

function Top
a = GetValue(1);
b = GetValue(2);
s = Sum2Num(a,b)
end

function val = GetValue(n)
val = n;
end


function s = Sum2Num(a,b)
s = a + b;
end


How to read a formatted text (ASCII) file in Matlab?

1. If a text (ASCII) file contains data with different data types then 'load' command is not useful. Then, one option can be to use TEXTSCAN with the format specifier as shown below.


fid = fopen(fileName)
formatSpec = '%f %s %f %f'

data = textscan(fid,formatSpec,3600*1,'Delimiter',',','CollectOutput', true);

fclose(fid)

The output will be a cell array. Each cell contains each of the elements. If you...
CollectOutput - will collect the consecutive data of same type in a single cell array.
Delimiter - lets you specify the delimiter used in the text file. If do not use the 'Delimiter' option then space is treated as the delimiter.


2. Many examples are given in the Matlab link given above:


Friday 21 November 2014

How to stop Matlab at the point of error for debugging?

1. There is an useful command 'dbstop'. Check help for more details.

2. To stop the control at the point of error in part of the function or the script use

dbstop if error

3. This has to be only once in the file. It gets applied to all the files.

4. Quite useful. Avoids rerunning the script to look for the reasons for the errors. It stops the control at the line where error is encountered and retains the work-space variables that belong to that function.


Thursday 13 November 2014

How to design a simple low pass filter (LPF) in Matlab?

%% LPF filter specifications

fSamp = 45e6;

Fpass = 1.25e6/(fSamp/2);          % Passband Frequency
Fstop = 2.5e6/(fSamp/2);        % Stopband Frequency
Apass = 1;           % Passband Ripple (dB)
Astop = 80;          % Stopband Attenuation (dB)
match = 'passband';  % Band to match exactly

% Construct an FDESIGN object and call its BUTTER method.
%d  = fdesign.lowpass(Fpass, Fstop, Apass, Astop, fSamp);
d=fdesign.lowpass('Fp,Fst,Ap,Ast',Fpass,Fstop,Apass,Astop);
designmethods(d)
Hd = design(d,'equiripple');
fvtool(Hd); %shows the response

%Using filter

inpSignalFiltered = filter(Hd, inpSignal);




Wednesday 6 August 2014

How to add colorbar ( adjust its size and label) to a plot in Matlab?

1. Use 'colorbar'

c=colorbar('vert'); % vert adds colorbar to z-axis values. 'horiz' adds color bar to x-values.

2. Reduce the size of the colorbar

xa1=get(gca,'position');
xa=get(c,'Position');
xa(3)=0.03;
xa(4)=0.65;

set(c,'Position',xa)

3. Adding title to the colorbar (this sits on top)

title(c,'Hz')

4. Adding ylabel to it (along vertical axis)

ylabel('Error')

Example:


Thursday 29 May 2014

How to get all the elements of a structure into single array, in Matlab?

 %Structure definition


% Case 1: Each element in the structure has same number of columns but different
% number of rows.

a(1).res = ones(4,2)*100;
a(2).res = ones(2,2)*2;
a(3).res = ones(3,2)*40;
a(4).res = ones(5,2)*253;

% For Case 1, [a.res] = horzcat will not work because it requires rows of
% all the elements in the structure to be equal. However, we can do
% vertical concatenation as columns of all the elements are equal.


vertCatAres = vertcat(a.res)

%Case 2: Each element in the structure has same number of rows but
%different number of columns

a(1).res = ones(2,4)*100;
a(2).res = ones(2,2)*2;
a(3).res = ones(2,3)*40;
a(4).res = ones(2,5)*253;
horzCatAres = [a.res] % horzcat(a.res)

%Results
vertCatAres =
   100   100
   100   100
   100   100
   100   100
     2     2
     2     2
    40    40
    40    40
    40    40
   253   253
   253   253
   253   253
   253   253
   253   253

horzCatAres =
   100   100   100   100     2     2    40    40    40   253   253   253   253   253
   100   100   100   100     2     2    40    40    40   253   253   253   253   253

Friday 28 March 2014

How to go to a particular byte in a file in Matlab?

1. Use fseek

2. Note: 
a. In file, byte numbering starts from 0. Therefore, if we want to go to Nth byte, then we have to move the pointer by (N-1) bytes.

fseek(fid, N-1,'bof'); %skips (N-1) bytes and points to Nth byte of the data.

b. In file, byte numbering starts from 0. Therefore, if we want to remove N bytes, move the pointer by N bytes.

fseek(fid, N,'bof'); % skips N bytes and points to N+1 byte of the data.

c. The N has to be an integer.

3. Example: 

status =  fseek(fid3, floor(numberOfBytesToSkipFromBOFFor20msIntAlignment), 'bof')

ferror(fid3) ; % ferror gives more information if there is any error during fseek operation.

status = 1 means successful operation else there is an error.

Tuesday 11 March 2014

How to put legends only to some of the lines in a Matlab Plot?

1. Use the 'plot' handles.

x1 = 1:10;
y1blue = 2 .* x1;


y2green = 2 .* y1;

y3red = y2 + 5;

h1 = plot(x1,y1blue ,'b');
hold on;

h2 = plot(y2green,'g');
h3 = plot(y3red,'r');

%Legend all the line plots
legend('Y1blue','Y2green', 'Y3red');


%Legend only Y1 and Y2
legend([h1 h2],{'Y1','Y2'});

2. Another method to control the 'legend content' is given HERE


3. How to use numbers in the legend?

legend({['PRN' num2str(1)],['PRN' num2str(2)]});

Monday 17 February 2014

Plot the spectrum of the real sinuoids and complex sinuoids.

Magnitude spectrum:
Magnitude spectrum of the sinusoids

Phase spectrum:

Phase spectrum (-fs/2 to +fs/2)

3D-Plots of Sinusoids: [note that Sin(-wt) should be read as Sin(wt)]
Sinusoids on Real and Imaginary axes

Complex Sinusoids 


%Code used to generate above plots

%% Clear the workspace
clear all;
close all;
clc;
home;

%% Parameters
fSamp = 1e3; % Hz
ts = 1/fSamp;
T = 1; % seconds - Data duration
t = 0:ts:T-ts;

fSig = 10; % Hz - Signal frequency
%% Real - sinusoids
% 1. Sine
s = sin(2*pi*fSig*t);

% 2. Co-sine
c = cos(2*pi*fSig*t);

% 3. Complex exponential
compExpoPositive = exp(1i*2*pi*fSig*t);
compExpoNegative = exp(-1i*2*pi*fSig*t);

%% Their Magnitude spectrum
df = fSamp / length(s);
freqIndex = -fSamp/2:df:fSamp/2 - df;

sFft = fftshift(fft(s)) ./ length(s);
% phaseSine = unwrap(angle((sFft))); %output is in radians but continuos (unwrap)
phaseSine = phase(sFft); %equal to unwrap and angle

cFft = fftshift(fft(c))./ length(s);
phaseCosine = unwrap(angle((cFft)));

fftCompPositive = fftshift(fft(compExpoPositive))./ length(s);
phaseCompPositive = phase(fftCompPositive);
fftCompNegative = fftshift(fft(compExpoNegative))./ length(s);
phaseCompNegative = phase(fftCompNegative);
% Sine
figure;
subplot(3,1,1);
stem(freqIndex,real(sFft),'b*'); hold on;
stem(freqIndex,imag(sFft),'ro');


xlim([-15 15]);
title('Sinusoid''s magnitude spectrum','color','b');
xlabel('Frequency (Hz)');
ylabel('Mag.');
legend('Real','Imag','location','best');

% Cosine
subplot(3,1,2);
stem(freqIndex,real(cFft),'g-'); hold on;
stem(freqIndex,imag(cFft),'m^');
xlim([-15 15]);
title('Co-sinusoid''s magnitude spectrum','color','b');
xlabel('Frequency (Hz)');
ylabel('Mag.');
legend('Real','Imag','location','best');

% Complex sinusoids
subplot(3,1,3);
stem(freqIndex,real(fftCompPositive),'b*'); hold on;
stem(freqIndex,imag(fftCompPositive),'ro');

stem(freqIndex,real(fftCompNegative),'g*');%,'markersize',5);
stem(freqIndex,imag(fftCompNegative),'mo');%,'markersize',5);
xlim([-15 15]);
title('Complex exponentials / Complex sinuoids','color','b');
xlabel('Frequency (Hz)');
ylabel('Mag. (arbitrary unit)');
legend('Real(exp(jwt))','Imag(exp(jwt))','Real(exp(-jwt))','Imag(exp(-jwt))','location','best');

% Setting plot properties
set(findall(gcf,'Type','line'),'linewidth',2);
set(findall(gcf,'type','axes'),'fontsize',12);
set(findall(gcf,'Type','text'),'FontSize',12);
saveas(gcf,'MagSpectrumFreqDomain.jpg');




%% Their Phase spectrum
figure;
plot(freqIndex,phaseSine*180/pi,'b');
hold on;
plot(freqIndex,phaseCosine*180/pi,'r');

plot(freqIndex,phaseCompPositive*180/pi,'g');
plot(freqIndex,phaseCompNegative*180/pi,'m');


xlim([-15 15]);
% Setting plot properties
title('Phase plot');
xlabel('Frequency (Hz)');
ylabel('Phase (degrees)');



legend('Sine','Co-sine','exp(+jwt)','exp(-jwt)','location','best');
set(findall(gcf,'Type','line'),'linewidth',2);
set(findall(gcf,'type','axes'),'fontsize',12);
set(findall(gcf,'Type','text'),'FontSize',12);

saveas(gcf,'PhaseSpecFreqDomain.jpg');




***************************************************************** 







%% 3D plot of the real and complex sinusoids

% Real sinusoids

% Complex sinusoids

% Code used to generate 3-D plots

%% Clear the workspace

clc;
clear all;
close all;

figure;

%% Plotting components of Co-Sinusoids

x=ones(1,6) .* 0.5;

y=0:0.1:0.5;

z=zeros(1,length(x));

plot3(x,y,z,'g');
hold on;

x=ones(1,6) .* -0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));

plot3(x,y,z,'g+');

%% Plotting components of Sinusoids

x=ones(1,6) .* 0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));

plot3(x,z,-y,'m');

x=ones(1,6) .* -0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));
plot3(x,z,y,'m+');

%% Plotting axis lines

xL = [-1 1];%xlim;
yL = xL;
zL = yL;
line([0 0], yL, [0 0],'linewidth',3,'color','y'); %x-axis
line(xL, [0 0], [0 0],'linewidth',3,'color','k'); %y-axis
line([0 0], [0 0],zL,'linewidth',3,'color','b'); %y-axis
hold on;
grid on;

%%
axis equal;
set(findall(gcf,'Type','line'),'linewidth',2);
legend('Cos(wt) @ w', 'Cos(wt) @ -w', 'Sin(wt) @ w','Sin(-wt) @ -w');
ylabel('Real');
xlabel('Frequency (rad)');
zlabel('Imaginary')
view(36,25);
saveas(gcf,'RealSinusoids.jpg');

%% Plotting complex exponentials

figure;

x=ones(1,6) .* +0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));
plot3(x,y,z,'r');
hold on;
x=ones(1,6) .* -0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));
plot3(x,y,z,'b');

%
x=ones(1,6) .* 0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));
plot3(x,z,-y,'r+');

x=ones(1,6) .* -0.5;
y=0:0.1:0.5;
z=zeros(1,length(x));
plot3(x,z,y,'b+');

%% Plotting lines
xL = [-1 1];%xlim;
yL = xL;
zL = yL;
line([0 0], yL, [0 0],'linewidth',3,'color','y'); %x-axis
line(xL, [0 0], [0 0],'linewidth',3,'color','k'); %y-axis
line([0 0], [0 0],zL,'linewidth',3,'color','b'); %y-axis
hold on;
grid on;
axis equal;
set(findall(gcf,'Type','line'),'linewidth',2);

legend('exp(+jwt) ','exp(-jwt)','-j*exp(+jwt)','j*exp(-jwt)');
ylabel('Real');
xlabel('Frequency (rad)');
zlabel('Imaginary')
whitebg('white');
% whitebg('blue')
view(36,25);
saveas(gcf,'ComplexExpo.jpg');    






Thursday 6 February 2014

Euler's formula and How to plot the XY axis at the center of the plot?


%% Clearing workspace
clc;
clear all;
close all;
home;

%% Just a circle -- Euler's circle
f = 1; tm = 1/f;
numberOfSamplesPerCycle = 6; % define the sides of the polygon, if large the polygon appears as circle. If = 6, we get hexagon.

dt = tm / numberOfSamplesPerCycle;

t = 0:dt:(1/f); % the number of sides in the circle depend on the 'dt' step - this define the sample time interval

s = exp(1i*2*pi*f*t);


%% Another way to plot a circle -- Trigonometric relation with complex numbers

% We can get the same polygon using this below code:
r = abs(s);
theta = phase(s);
x = r .* cos(theta);
y = r .* sin(theta);
figure; plot(s,'g'); hold on; plot(x,y,'r*');% hold on;
legend('Euler''s style', 'polar style','Location','best');
xlabel('Real part');
ylabel('Imaginary part');
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 12, 'fontWeight', 'normal')
set(findall(gcf,'Type','line'),'linewidth',2);
saveas(gcf,'fig1.jpg');
Circle plot using Euler's formula and polar form formula

%% To draw the XY axes at the centre of the plot
xL = xlim;
yL = ylim;
line([0 0], yL); %x-axis
line(xL, [0 0]); %y-axis


xlabel('Real part');
ylabel('Imaginary part');
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 12, 'fontWeight', 'normal')
set(findall(gcf,'Type','line'),'linewidth',2);
saveas(gcf,'fig2.jpg');

Plotting Center XY coordinate lines

%% To draw a circle
f = 1; tm = 1/f;
numberOfSamplesPerCycle = 60; % define the sides of the polygon, if large the polygon appears as circle. If = 6, we get hexagon.

dt = tm / numberOfSamplesPerCycle;

t = 0:dt:(1/f); % the number of sides in the circle depend on the 'dt' step - this define the sample time interval
s = exp(1i*2*pi*f*t);

figure; plot(s,'r');

xL = xlim;
yL = ylim;
line([0 0], yL); %x-axis
line(xL, [0 0]); %y-axis

xlabel('Real part');
ylabel('Imaginary part');
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 12, 'fontWeight', 'normal')
set(findall(gcf,'Type','line'),'linewidth',2);

saveas(gcf,'fig3.jpg');
Smooth Cirle



%% Concentric circles
r1 = 0:0.1:50;
theta1 = r1;
x1 = r1 .* cos(theta1);
y1 = r1 .* sin(theta1);

figure;plot(x1,y1,'g')
xL = xlim;
yL = ylim;
line([0 0], yL); %x-axis
line(xL, [0 0]); %y-axis

xlabel('Real part');
ylabel('Imaginary part');
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 12, 'fontWeight', 'normal')
set(findall(gcf,'Type','line'),'linewidth',2);
saveas(gcf,'fig4.jpg');
Spiral plot using Euler's formula


%% Complex plane: s-plane
sigma = 0:0.1:1;
jOmega = sigma;
s = sigma + 1i * jOmega;
figure;plot(s)
xlabel('Real part');
ylabel('Imaginary part');
set(findall(gcf, '-property', 'FontSize'), 'FontSize', 12, 'fontWeight', 'normal')
set(findall(gcf,'Type','line'),'linewidth',2);

saveas(gcf,'fig5.jpg');
S-plane plot of a complex number