% SURAJ BASTOLA %
% ==========CHECKING WHETHER SYSTEM IS STABLE OR NOT ==================
% H(z) = (2z+0.5)/(z^2+0.5z-1)
b = [2 0.5];
a=[1 0.5 -1];
sys = tf(b,a);
p = pole(sys);
z = zero(sys);
fprintf("Poles are = %d\n",p);
if(p<=1)
fprintf("System is stable\n");
else
fprintf("System is unstable\n");
end
fprintf("Zeros are = %d\n",z);
subplot(2,1,1)
impz(b,a)
subplot(2,1,2)
zplane(b,a)
%======= SURAJ BASTOLA ========================%
%========== Assuming two obstacles with one obstacle loss equal to 1/3 of another Obstacle Loss ==========
clc;
u=-4:0.001:4;
for (i=1:length(u))
if (u(i)<=-1)
f1(i) = 0;
elseif(u(i)<=0)
f1(i) = 0.5-0.62*u(i);
elseif(u(i)<=1)
f1(i) = 0.5*exp(-0.95*u(i));
elseif(u(i)<=2.4)
f1(i) = 0.4-(sqrt(0.1184-(0.38-0.1*u(i))^2));
else
f1(i) = (0.225079)./u(i);
end
end
% =========== Single Knife-edge Diffraction Loss==============
L1 = 20*log10(abs(f1));
disp(L1);
subplot(2,1,1)
plot(u,L1,'g')
grid on
title('Single Knife-Edge Diffraction')
% ================= Multiple Knife Edge Diffraction with Epstein-Peterson Method======================
nobstac = 2;
Lsum = L1+(L1/3);
L2 = (Lsum)+(nobstac*3.9);
disp(L2);
subplot(2,1,2)
plot(u,L2,'r')
grid on
title('Multiple Knife Edge Diffraction')
xlabel('Fresnel Diffraction Parameter')
ylabel('Diffraction Loss (dB)')
% =================== SURAJ BASTOLA% ===========
% ============== Calculating Diffraction Loss ==============
clc;
clear all;
close all;
d1 = 1000;
d2 = 3480;
d3 = 1830;
d4 = 2650;
d5 = 3280;
d6 = 1200;
c = 3E8;
f = 8E8;
l = c/f;
fprintf('Wavelength = %d\n',l);
h1 = sqrt(l*d1*d2/(d1+d2));
h2 = sqrt(l*d3*d4/(d3+d4));
h3 = sqrt(l*d5*d6/(d5+d6));
fprintf('Height of knife edge above LOS of knife edge = %d %d %d \n',h1,h2,h3);
% =========== Obstacle 1 ==================
v1 = h1*(sqrt(2*(d1+d2)/(l*d1*d2)));
fprintf('First Fresnel-Kirchoff value = %d\n',v1);
if (v1<=-1)
f1 = 0;
elseif(v1<=0)
f1 = 0.5-0.62*v1;
elseif(v1<=1)
f1 = 0.5*exp(-0.95*v1);
elseif(v1<=2.4)
f1 = 0.4-(sqrt(0.1184-(0.38-0.1*v1)^2)); % It is our desired equation as v = 1.414214
else
f1 = (0.225079)./v1;
end
fprintf('Diffraction Loss =%d/n',f1); % Diffraction Loss
y1 = (abs(f1))^2;
G1 = 10* log10(y1); % Diffraction Loss in dB
fprintf('Total Diffraction Loss for X = %d\n',G1);
%======================Obstacle 2 =====================================
v2 = h2*(sqrt(2*(d3+d4)/(l*d3*d4)));
fprintf('First Fresnel-Kirchoff value = %d\n',v2);
%%%%%% Calculating Diffraction LOss in dB %%%%%%%%%%
if (v2<=-1)
f2 = 0;
elseif(v2<=0)
f2 = 0.5-0.62*v2;
elseif(v2<=1)
f2 = 0.5*exp(-0.95*v2);
elseif(v2<=2.4)
f2 = 0.4-(sqrt(0.1184-(0.38-0.1*v2)^2)); % It is our desired equation as v = 1.414214
else
f2 = (0.225079)./v2;
end
fprintf('Diffraction Loss =%d/n',f1); % Diffraction Loss
y2 = (abs(f2))^2;
G2 = 10* log10(y2); % Diffraction Loss in dB
fprintf('Total Diffraction Loss for X = %d\n',G2);
%======================Obstacle 3 =====================================
v3 = h3*(sqrt(2*(d5+d6)/(l*d5*d6)));
fprintf('First Fresnel-Kirchoff value = %d\n',v3);
%%%%%% Calculating Diffraction LOss in dB %%%%%%%%%%
if (v3<=-1)
f3 = 0;
elseif(v3<=0)
f3 = 0.5-0.62*v3;
elseif(v3<=1)
f3 = 0.5*exp(-0.95*v3);
elseif(v3<=2.4)
f3 = 0.4-(sqrt(0.1184-(0.38-0.1*v3)^2)); % It is our desired equation as v = 1.414214
else
f3 = (0.225079)./v3;
end
fprintf('Diffraction Loss =%d/n',f3); % Diffraction Loss
y3 = (abs(f3))^2;
G3 = 10* log10(y3); % Diffraction Loss in dB
fprintf('Total Diffraction Loss for X = %d\n',G3);
%======================Epstein Peterson Method with Foose Correction=====================================
nobst = 3;
L = G1+G2+G3+(nobst*3.91);
fprintf("Total diffraction Loss is %d\n", L);
% SURAJ BASTOLA %%%%%%%%%%%%
%%% RECORDING AND SAVING RECORDED AUDIO %%%%%%
% =======================Recording Audio and Saving it to some location================
audio = audiorecorder(8000,8,1); % y = audiorecorder(Fs,nbits,channels)
fprintf('start speaking\n')
%%%% Stopping record after 10 seconds %%%%%%%%
recordblocking(audio,10); % recordblocking(recorderObj, length)
fprintf('End of recording\n')
%%%%% Getting data from audiorecorder Object and saving it to specified folder location %%%%%%%%%%%
myrec= getaudiodata(audio); % y = getaudiodata(recorder)
audiowrite('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav', myrec, 8000);
fprintf('successfully saved the recorded audio');
% ========================Loading recorded audio by retrieving from saved location ===================
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav');
sound(y,8000); % Playing retrieved audio
fprintf('successfully played the retrieved recorded audio\n')
% ===== Plotting the audio %%%%%%%%%%%
plot(y)
title('Recorded Audio of 10 seconds')
% %%%%% Orthogonality %%%%%%%%%%
%%%%%%% SURAJ BASTOLA %%%%%%%%%%%%%%
N = 4;
df = 2;
f = 100:df:20+(N-1)*df;
T = 1/df; % different values of symbol duration
t = linspace(0,T,1000);
x = zeros(N, length(t));
for i = 1:length(f)
x(i,:) = sin(2*pi*f(i)*t);
subplot(N*2,1,i);
plot(t,x(i,:));
end
for i = 1:N
for j = i+1:N
if i==1
subplot(N*2,1,N+j);
plot(t,x(i,:).*x(j,:));
end
fprintf("Correlation(f%d,f%d) = %f \n",i,j,mean(x(i,:)),x(j,:));
end
end
%%%%%%%%%%% SURAJ BASTOLA %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% https://lnkd.in/gvdzPWX %%%%%%%%%%%%%
%%%%%%%%% RICEAN AND RAYLEIGH SIGNAL FADING %%%%%%%%%
clc
close all;
%N = 10; % number of samples (bits) for transmission
k=1; %Ricean factor
SNR = 8; % SNR in dB
SNR1 = 10^(SNR./10); % or, SNR1 = 10^(SNR/10)
%disp(SNR1);
% Signal bit creation
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB.wav');
subplot(2,2,1)
plot(y)
title('Original Audio Signal')
Sig_Bit = dec2bin( typecast( single(y), 'uint8'), 8 ) - '0';
Sig_Bit = reshape(Sig_Bit,1,[]);
N = 512;
SignalX = Sig_Bit(20000:20000+N-1);
%Sig_Bit = num2bin(y);
%Sig_Bit = randn(1,N)>0.5;
fprintf('Signal Created');
%disp(Sig_Bit);
%BPSK symbol
raw_signal = 2*SignalX-1;
fprintf('BPSK SYMBOL');
subplot(2,2,2)
plot(raw_signal)
title('Raw BPSK Signal')
%disp(raw_signal);
n = 1./sqrt(2)*(randn(1,N) + 1j*randn(1,N));%Noise creation(AWGN, 0dB variance )
%================== RAYLEIGH CHANNEL ================================
ch_ray = 1./sqrt(2)*(randn(1,N)+1i*randn(1,N));
% Faded signal with Rayleigh channel
Sig_fad_ray = ch_ray.* raw_signal + SNR1.*n;
subplot(2,2,3)
plot(Sig_fad_ray)
title('Faded Rayleigh Signal')
%====================================================
% equalization (ideal)
Sig_fad_ray1 = Sig_fad_ray./ch_ray;
fprintf('Rayleigh Faded Signal Received at the receiver');
%disp(Sig_fad_ray1);
fprintf('Absolute value of Faded Signal Received at the receiver');
%disp(abs(Sig_fad_ray1));
%================== RICEAN CHANNEL ================================
% Ricean Channel
ch_ric = sqrt(k/(k+1))+sqrt(1/(k+1))*(1/sqrt(2))*(randn(1,N)+1j*randn(1,N));
Sign_fad_ric = ch_ric.*raw_signal+10^(-SNR/20)*n;
% equalization (ideal)
Sig_fad_ric1 = Sign_fad_ric./ch_ric;
fprintf('Faded Signal Received at the receiver');
%disp(Sig_fad_ric1);
fprintf('Absolute value of Faded Signal Received at the receiver');
%disp(abs(Sig_fad_ric1));
subplot(2,2,4)
plot(Sig_fad_ric1)
title('Faded Ricean Signal')
% receiver - hard decision decoding
for i = 1: N
if real(Sig_fad_ray1(i)) > 0
Sig_ray_Rayleigh(i) = 1;
else
Sig_ray_Rayleigh(i) = 0;
end
end
% receiver - hard decision decoding
for i = 1: N
if real(Sig_fad_ric1(i)) > 0
Sig_ray_Rec(i) = 1;
else
Sig_ray_Rec(i) = 0;
end
end
% Measuring the errors: Compare original signal with
Err_num_ray1 = size(find([SignalX- Sig_ray_Rayleigh]),2);
BER_ray1 = Err_num_ray1 /N;
fprintf('Rayleigh BER = ');
disp(BER_ray1);
Recian = size(find([SignalX- Sig_ray_Rec]),2);
BER_Recian = Recian /N;
fprintf('Ricean BER= ');
disp(BER_Recian);
%%%%%%%%%% SURAJ BASTOLA %%%%%%%%%%
M = 8; % for 8 BPSK
N = 64;% this means we are using 64 point FFT/IFFT
data = randi(M,1,N)-1;
subplot(411);
myplot(stem(data),'Data',3);
X = pskmod(data,M);
figure(2);
subplot(111);
myplot(plot(X,'o'),'8PSK',2);
% IFFT with real result for Baseband
% for broad band go directly to line 17 from 10 with code x = ifft(X, N);
X1 = X;
X2 = conj(X(N:-1:2));
X1(1) = real(X(1));
X = [X1, imag(X(1)), X2];
x = ifft(X, N*2);
figure(1)
subplot(412);
myplot(plot(1:128,x),'IFFT', 1);
% add circular prefix
delay = 10;
xt = [x(2*N-delay+1:2*N),x];
subplot(413);
myplot(plot(1:(2*N+delay),xt),'CP',2);
hold on
myplot(plot(1+delay:(2*N+delay),x),'OFDM',1);
% remove Cyclic Prefix
r = xt((1+delay):(2*N+delay));
%fft
r2 = fft(r,N*2);
% undo complex number modification
r3 = r2(1:N);
r3(1) = r2(1)+1i*r2(N+1);
data2 = pskdemod(r3,M);
subplot(414);
myplot(stem(data2),'Received',4);
disp([data;data2]);
%======= SURAJ BASTOLA ========================%
%==========Plotting
Frenel Diffraction Paramter with Diffraction Loss%=====
clc;
clear all;
u=-5:0.05:5;
for (i=1:length(u))
if (u(i)<=-1)
f1(i) = 0;
elseif(u(i)<=0)
f1(i) = 0.5-0.62*u(i);
elseif(u(i)<=1)
f1(i) = 0.5*exp(-0.95*u(i));
elseif(u(i)<=2.4)
f1(i) = 0.4-(sqrt(0.1184-(0.38-0.1*u(i))^2));
else
f1(i) = (0.225079)./u(i);
end
end
G = 20*log(abs(f1));
plot(u,G,'r')
xlabel('Fresnel Diffraction Parameter')
ylabel('Diffraction Loss (dB)')
title('Frenel Diffraction Paramter Analysis with Diffraction Loss')
% =================== SURAJ BASTOLA% ===========
%==================THANKS TO MATLAB =========================
% ============== Calculating Diffraction Loss ==============
clc;
clear all;
close all;
d1 = 20000;
d2 = 36500;
c = 3E8;
f = 9E8;
l = c/f;
fprintf('Wavelength = %d\n',l);
h = sqrt(l*d1*d2/(d1+d2));
fprintf('Height of knife edge above LOS of knife edge = %d\n',h);
%%%%%%%% Calculating first Fresnel-Kirchhoff value %%%%%%%
v = h*(sqrt(2*(d1+d2)/(l*d1*d2)));
fprintf('First Fresnel-Kirchoff value = %d\n',v);
%%%%%% Calculating Diffraction LOss in dB %%%%%%%%%%
if (v<=-1)
f = 0;
elseif(v<=0)
f = 0.5-0.62*v;
elseif(v<=1)
f = 0.5*exp(-0.95*v);
elseif(v<=2.4)
f = 0.4-(sqrt(0.1184-(0.38-0.1*v)^2)); % It is our desired equation as v = 1.414214
else
f = (0.225079)./v;
end
fprintf('Diffraction Loss =%d/n',f); % Diffraction Loss
y = (abs(f))^2;
G = 10* log10(y); % Diffraction Loss in dB
fprintf('Total Diffraction Loss in dB = %d\n',G);
======================================================
Output:
Wavelength = 3.333333e-01
Height of knife edge above LOS of knife edge = 6.562610e+01
First Fresnel-Kirchoff value = 1.414214e+00
Diffraction Loss =1.520479e-01/nTotal Diffraction Loss in dB = -1.636039e+01
clc;
% -------- SURAJ BASTOLA ------------%
% -------- Refractivity Calculation-------------%
h0 = 7.350; % height above sea-level
N0 = 315;% Refractivity above sea-level
h = 0:0.5:8.848; % Mt. Everest = 8.848km i.e. height above earth surface
n = 1+ (N0*10^-6)*exp(-h/h0); % refractivity above earth surface
display(h);
display(n);
plot(h,n) % plotting refractivity vs height above earth surface
grid on;
xlabel("Height above earth surface");
ylabel("Refractivity above earth surface");
title("REFRACTIVITY CALCULATION FROM FOOT TO TOP OF MT.Everest");
%%%%%%%%% QAM Analysis %%%%%%%%%%%%
clear all;
clc
close all;
M = 256;
k = 0:M-1;
Fs = 8000;
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav');
data = dec2bin( typecast( single(y), 'uint8'), 8 ) - '0';
A = qammod(data,M);% QAM symbol generation
x = zeros(1,M);
for i = k
sum = 0;
for m = k
sum = sum + A(m+1)*exp(2i*pi*m*i/M);
end
x(i+1) = sum;
end
x2 = ifft(A,M);
subplot(211);
%plot(abs(x2))
hist(real(x2),5);
title('Real Part')
subplot(212)
hist(imag(x2),5);
%plot(imag(x2))
title('Imaginary Part')
sgtitle('QAM signal generation')
% SURAJ BASTOLA %%%%%%%%%%%%
%%% Low Pass Filter %%%%%%%%%%
% IN MATLAB, Fs takes default value of 8000 %%%%%%%%%%%
clc;
clear all;
close all;
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav');
z = abs(fft(y));
disp(z);
figure(1)
subplot(2,3,1)
plot(y)
title('Original Recorded Audio')
t = (0:length(y)-1)/Fs;
% ================ Low pass filter using Chebyshev Window==============
% I have used 34th order FIR low pass filter
% which has cut-off frequency0.48 and Chebyshev window with 30dB of ripple
blo = fir1(34,0.48,chebwin(35,30));
outlo = filter(blo,1,y);
ys = ylim;
subplot(2,3,2)
plot(t,outlo)
title('Lowpass Filtered Signal')
xlabel('Time (s)')
ylim(ys)
% ================ Low pass filter using Hamming Window==============
ord = 46;
low = 0.4;
bnd = [0.6 0.9];
bM = fir1(ord, [low bnd]);
Hamming = filter(bM,1,y);
subplot(2,3,3)
plot(t, Hamming)
title('Hamming window')
%=================== Low pass filter using Hamming Window==========
hM = fir1(ord,[low bnd],'DC-0',hann(ord+1));
Hanning = filter(hM,1,y);
subplot(2,3,4)
plot(t, Hanning)
title('Hanning window')
hfvt = fvtool(blo,1,bM,1,hM,1);
hfvt = fvtool(outlo,1,Hamming,1,Hanning,1);
legend(hfvt,'Chebwin','Hamming','Hanning')
% SURAJ BASTOLA %%%%%%%%%%%%
%%% COMPRESSING THE RECORDED AUDIO USING A LAW%%%%%%
% =============== Syntax: out = compand(in,A,v,'A/compressor') ==========================================
% IN MATLAB, Fs takes default value of 8000
clc;
close all;
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav');
z = abs(fft(y));
disp(z);
figure(1)
subplot(2,2,1)
plot(y)
title('Recorded Audio')
subplot(2,2,2)
plot(z)
title('FFT of Recorded Audio')
A = 87.6
compressed = compand(y,A,max(y),'A/compressor');
x = abs(fft(compressed));
disp(x);
subplot(2,2,3)
plot(compressed)
title('Compressed Audio')
subplot(2,2,4)
plot(x)
title('FFT of Compressed Audio')
sgtitle('Compressing audio using A-law');
% SURAJ BASTOLA %%%%%%%%%%%%
%%% COMPRESSING THE RECORDED AUDIO USING Mu-LAW%%%%%%
% IN MATLAB, Fs takes default value of 8000
% ===================== Syntax: out = compand(in,Mu,v,'mu/compressor') =======================================================
clc;
close all;
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB COding\audioRecorded.wav');
z = abs(fft(y));
disp(z);
figure(1)
subplot(2,2,1)
plot(y)
title('Recorded Audio')
subplot(2,2,2)
plot(z)
title('FFT of Recorded Audio')
% Mu = 255
compressed = compand(y,255,max(y),'mu/compressor'); % out = compand(in,Mu,v) implements a µ-law compressor for the input vector in. Mu specifies µ, and v is the input signal's maximum magnitude. out has the same dimensions and maximum magnitude as in.
x = abs(fft(compressed));
disp(x);
subplot(2,2,3)
plot(compressed)
title('Compressed Audio')
subplot(2,2,4)
plot(x)
title('FFT of Compressed Audio')
sgtitle('Compressing audio using Mu- law');
%%%%%%%%%%% SURAJ BASTOLA %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% https://lnkd.in/gvdzPWX %%%%%%%%%%%%%
%%%%%%%%%% MIT183249 %%%%%%%%%%%%%%
clear all;
clc
close all;
%N = 10; % number of samples (bits) for transmission
k=1; %Ricean factor
SNR = 8; % SNR in dB
SNR1 = 10^(SNR./10); % or, SNR1 = 10^(SNR/10)
%disp(SNR1);
% Signal bit creation
[y,Fs] = audioread('C:\Users\suraj\OneDrive\Documents\MIT COLLEGE\Software Defined Radio Communication\MATLAB.wav');
Sig_Bit = dec2bin( typecast( single(y), 'uint8'), 8 ) - '0';
Sig_Bit = reshape(Sig_Bit,1,[]);
N = 512;
SignalX = Sig_Bit(20000:20000+N-1);
%Sig_Bit = num2bin(y);
%Sig_Bit = randn(1,N)>0.5;
fprintf('Signal Created');
%disp(Sig_Bit);
%BPSK symbol
raw_signal = 2*SignalX-1;
fprintf('BPSK SYMBOL');
%disp(raw_signal);
n = 1./sqrt(2)*(randn(1,N) + 1j*randn(1,N));%Noise creation(AWGN, 0dB variance )
%================== RAYLEIGH CHANNEL ================================
ch_ray = 1./sqrt(2)*(randn(1,N)+1i*randn(1,N));
% Faded signal with Rayleigh channel
Sig_fad_ray = ch_ray.* raw_signal + SNR1.*n;
%====================================================
% equalization (ideal)
Sig_fad_ray1 = Sig_fad_ray./ch_ray;
fprintf('Rayleigh Faded Signal Received at the receiver');
%disp(Sig_fad_ray1);
fprintf('Absolute value of Faded Signal Received at the receiver');
%disp(abs(Sig_fad_ray1));
%================== RICEAN CHANNEL ================================
% Ricean Channel
ch_ric = sqrt(k/(k+1))+sqrt(1/(k+1))*(1/sqrt(2))*(randn(1,N)+1j*randn(1,N));
Sign_fad_ric = ch_ric.*raw_signal+10^(-SNR/20)*n;
% equalization (ideal)
Sig_fad_ric1 = Sign_fad_ric./ch_ric;
fprintf('Faded Signal Received at the receiver');
%disp(Sig_fad_ric1);
fprintf('Absolute value of Faded Signal Received at the receiver');
%disp(abs(Sig_fad_ric1));
% receiver - hard decision decoding
for i = 1: N
if real(Sig_fad_ray1(i)) > 0
Sig_ray_Rayleigh(i) = 1;
else
Sig_ray_Rayleigh(i) = 0;
end
end
% receiver - hard decision decoding
for i = 1: N
if real(Sig_fad_ric1(i)) > 0
Sig_ray_Rec(i) = 1;
else
Sig_ray_Rec(i) = 0;
end
end
% Measuring the errors: Compare original signal with
Err_num_ray1 = size(find([SignalX- Sig_ray_Rayleigh]),2);
BER_ray1 = Err_num_ray1 /N;
fprintf('Rayleigh BER = ');
disp(BER_ray1);
Recian = size(find([SignalX- Sig_ray_Rec]),2);
BER_Recian = Recian /N;
fprintf('Ricean BER= ');
disp(BER_Recian);
% SURAJ BASTOLA % % ==========CHECKING WHETHER SYSTEM IS STABLE OR NOT ================== % H(z) = (2z+0.5)/(z^2+0.5z-1) b = [2 0.5]; a=[1...
