From ff1fd88553ad53b223f51518b534497db4274124 Mon Sep 17 00:00:00 2001 From: Ahxmedxx Date: Wed, 20 Dec 2023 12:09:25 +0530 Subject: [PATCH] Commiting all the Matlab Codes --- BW_IIR filter.m | 64 ++++++++++++++++++++++++++++++++++++ BW_IIR_filter.m | 64 ++++++++++++++++++++++++++++++++++++ Chebyshev_IIR.m | 33 +++++++++++++++++++ DFT.m | 26 +++++++++++++++ DFT2.m | 26 +++++++++++++++ FFT.m | 67 ++++++++++++++++++++++++++++++++++++++ FIR.m | 58 +++++++++++++++++++++++++++++++++ FIR1.m | 49 ++++++++++++++++++++++++++++ IDFT.m | 25 +++++++++++++++ IFFT.m | 73 ++++++++++++++++++++++++++++++++++++++++++ IIIR.m | 14 ++++++++ IIR.m | 15 +++++++++ IIR_basic.m | 15 +++++++++ SP_OPEN_ENDED.m | 73 ++++++++++++++++++++++++++++++++++++++++++ SSPP.m | 72 +++++++++++++++++++++++++++++++++++++++++ circular_convolution.m | 23 +++++++++++++ comb_func.m | 57 +++++++++++++++++++++++++++++++++ linear_circular.m | 14 ++++++++ linear_conv_fft.m | 53 ++++++++++++++++++++++++++++++ linear_convolution.m | 41 ++++++++++++++++++++++++ sampling_func.m | 49 ++++++++++++++++++++++++++++ signals.m | 54 +++++++++++++++++++++++++++++++ sine_func.m | 15 +++++++++ 23 files changed, 980 insertions(+) create mode 100644 BW_IIR filter.m create mode 100644 BW_IIR_filter.m create mode 100644 Chebyshev_IIR.m create mode 100644 DFT.m create mode 100644 DFT2.m create mode 100644 FFT.m create mode 100644 FIR.m create mode 100644 FIR1.m create mode 100644 IDFT.m create mode 100644 IFFT.m create mode 100644 IIIR.m create mode 100644 IIR.m create mode 100644 IIR_basic.m create mode 100644 SP_OPEN_ENDED.m create mode 100644 SSPP.m create mode 100644 circular_convolution.m create mode 100644 comb_func.m create mode 100644 linear_circular.m create mode 100644 linear_conv_fft.m create mode 100644 linear_convolution.m create mode 100644 sampling_func.m create mode 100644 signals.m create mode 100644 sine_func.m diff --git a/BW_IIR filter.m b/BW_IIR filter.m new file mode 100644 index 0000000..fe0f17d --- /dev/null +++ b/BW_IIR filter.m @@ -0,0 +1,64 @@ +%IIR BUTTERWORTH FILTER +syms s +kp=-2; +ks=-10; +wp=20; +ws=30; +j=sqrt(-1); +N=log((10^(-kp/10)-1)/(10^(-ks/10)-1))/(2*log(wp/ws)) +N=round(N)+1 +wc=wp/(10^(-kp/10)-1)^(1/(2*N)) +H1=1/(s+1); +H2=1/(s^2+1.414*s+1); +H3=1/((s^2+s+1)*(s+1)); +H4=1/((s^2+0.765368*s+1)*(s^2+1.84776*s+1)); +H5=1/((s+1)*(s^2+0.6180*s+1)*(s^2+1.6180*s+1)); + +if N==1 +H(s)=H1; +elseif N==2 +H(s)=H2; +elseif N==3 +H(s)=H3; +elseif N==4 +H(s)=H4; +elseif N==5 +H(s)=H5; +end +H(s)=H(s/wc) +H(s)=H(j*wp) +n=abs(H(s)) + +m=20*log(H(s)) +%IIR BUTTERWORTH FILTER +syms s +kp=-2; +ks=-10; +wp=20; +ws=30; +j=sqrt(-1); +N=log((10^(-kp/10)-1)/(10^(-ks/10)-1))/(2*log(wp/ws)) +N=round(N)+1 +wc=wp/(10^(-kp/10)-1)^(1/(2*N)) +H1=1/(s+1); +H2=1/(s^2+1.414*s+1); +H3=1/((s^2+s+1)*(s+1)); +H4=1/((s^2+0.765368*s+1)*(s^2+1.84776*s+1)); +H5=1/((s+1)*(s^2+0.6180*s+1)*(s^2+1.6180*s+1)); + +if N==1 +H(s)=H1; +elseif N==2 +H(s)=H2; +elseif N==3 +H(s)=H3; +elseif N==4 +H(s)=H4; +elseif N==5 +H(s)=H5; +end +H(s)=H(s/wc) +H(s)=H(j*wp) +n=abs(H(s)) + +m=20*log(H(s)) diff --git a/BW_IIR_filter.m b/BW_IIR_filter.m new file mode 100644 index 0000000..fe0f17d --- /dev/null +++ b/BW_IIR_filter.m @@ -0,0 +1,64 @@ +%IIR BUTTERWORTH FILTER +syms s +kp=-2; +ks=-10; +wp=20; +ws=30; +j=sqrt(-1); +N=log((10^(-kp/10)-1)/(10^(-ks/10)-1))/(2*log(wp/ws)) +N=round(N)+1 +wc=wp/(10^(-kp/10)-1)^(1/(2*N)) +H1=1/(s+1); +H2=1/(s^2+1.414*s+1); +H3=1/((s^2+s+1)*(s+1)); +H4=1/((s^2+0.765368*s+1)*(s^2+1.84776*s+1)); +H5=1/((s+1)*(s^2+0.6180*s+1)*(s^2+1.6180*s+1)); + +if N==1 +H(s)=H1; +elseif N==2 +H(s)=H2; +elseif N==3 +H(s)=H3; +elseif N==4 +H(s)=H4; +elseif N==5 +H(s)=H5; +end +H(s)=H(s/wc) +H(s)=H(j*wp) +n=abs(H(s)) + +m=20*log(H(s)) +%IIR BUTTERWORTH FILTER +syms s +kp=-2; +ks=-10; +wp=20; +ws=30; +j=sqrt(-1); +N=log((10^(-kp/10)-1)/(10^(-ks/10)-1))/(2*log(wp/ws)) +N=round(N)+1 +wc=wp/(10^(-kp/10)-1)^(1/(2*N)) +H1=1/(s+1); +H2=1/(s^2+1.414*s+1); +H3=1/((s^2+s+1)*(s+1)); +H4=1/((s^2+0.765368*s+1)*(s^2+1.84776*s+1)); +H5=1/((s+1)*(s^2+0.6180*s+1)*(s^2+1.6180*s+1)); + +if N==1 +H(s)=H1; +elseif N==2 +H(s)=H2; +elseif N==3 +H(s)=H3; +elseif N==4 +H(s)=H4; +elseif N==5 +H(s)=H5; +end +H(s)=H(s/wc) +H(s)=H(j*wp) +n=abs(H(s)) + +m=20*log(H(s)) diff --git a/Chebyshev_IIR.m b/Chebyshev_IIR.m new file mode 100644 index 0000000..4ddbb18 --- /dev/null +++ b/Chebyshev_IIR.m @@ -0,0 +1,33 @@ +%Chebyshev Filter +syms s +kp=-2; +ks=-20; +wp=1; +ws=1.3; +j=sqrt(-1); +e=sqrt(((1/(10^(kp/20)))^2)-1); +dp=1-(1/sqrt(1+(e^2))); +ds=10^(ks/20); +%discrimination factor +d=sqrt(((1-dp)^(-2)-1)/((ds^(-2))-1)); +%selectivity factor +k=wp/ws; +N=(acosh(1/d))/(acosh(1/k)); +N=round(N)+1; +%to find a and b + +a1=(1/2)*(((1+sqrt(1+e^2))/e)^(1/N)); +a2=(1/2)*(((1+sqrt(1+e^2))/e)^(-1/N)); +a=a1-a2; +b=a1+a2; + +for k=1:1:N +S(k+1)=(-a)*sin(((2*k)-1)*(pi/(2*N))); +end + +for k=1:1:N +w(k+1)=b*cos(((2*k)-1)*(pi/(2*N))); +end +for m=1:1:N +l(m+1)=(s-(S(m)+(j*w(m)))); +end \ No newline at end of file diff --git a/DFT.m b/DFT.m new file mode 100644 index 0000000..daab82d --- /dev/null +++ b/DFT.m @@ -0,0 +1,26 @@ +x = [1,2,3,4]; +N = length(x); +% k = (1:N) +X = zeros(1, N); +j = sqrt(-1); +disp("DFT of x(n) ") +%DFT of x(n) +for k = 0:N-1 + for n= 0:N-1 + X(k+1) = X(k+1) + x(n+1)*(exp(-j * n * k *2 *pi / N)) + end +end +ny = (1:N); +X_mag = abs(X); +X_phase = angle(X); +figure("Name","DFT") +subplot(2,1,1) +stem(ny,X_mag,"LineWidth",3) +xlabel("Frequency") +ylabel("Mag") +title("Magnitude plot") +subplot(2,1,2) +stem(ny, X_phase,"LineWidth",3) +xlabel("Frequency") +ylabel("Phase") +title("Phase plot") diff --git a/DFT2.m b/DFT2.m new file mode 100644 index 0000000..17d64fc --- /dev/null +++ b/DFT2.m @@ -0,0 +1,26 @@ +x = [1,2,3,4,5]; +N = length(x); +% k = (1:N) +X = zeros(1, N); +j = sqrt(-1); +disp("DFT of x(n) ") +%DFT of x(n) +for k = 0:N-1 + for n= 0:N-1 + X(k+1) = X(k+1) + x(n+1)*(exp(-j * n * k *(2 *pi / N))) + end +end +ny = (1:N); +X_mag = abs(X); +X_phase = angle(X); +figure("Name","DFT") +subplot(2,1,1) +stem(ny,X_mag,"LineWidth",3) +xlabel("Frequency") +ylabel("Mag") +title("Magnitude plot") +subplot(2,1,2) +stem(ny, X_phase,"LineWidth",3) +xlabel("Frequency") +ylabel("Phase") +title("Phase plot") \ No newline at end of file diff --git a/FFT.m b/FFT.m new file mode 100644 index 0000000..c04f9d5 --- /dev/null +++ b/FFT.m @@ -0,0 +1,67 @@ +syms Wnkm(k,m) Wnk(k) +x = [1 2 3 4 5 6 7 8]; +N = length(x); +f1 = (zeros(1,round(N/2))); +f2 = (zeros(1,round(N/2))); +y=1; +%Decimating into EVEN and ODD +for i=1:N + if rem(i,2)==0 + f2(y) = x(i); + y=y+1; + else + f1(y) = x(i); +end +end +j=sqrt(-1); +Wnkm(k,m) = exp((-j*4*pi*k*m)/N); +Wnk(k) = exp((-j*2*pi*k)/N); +F1 = (zeros(1,round(N/2))); +F2 = (zeros(1,round(N/2))); +%F1(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F1(k) = F1(k) + f1(m)*Wnkm(k-1,m-1); + end +end +%F2(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F2(k) = F2(k) + f2(m)*Wnkm(k-1,m-1); + end +end +F1; +F2; +X1 = (zeros(1,round(N/2))); +X2 = (zeros(1,round(N/2))); +%X1(k) = X(k) and X2(k) = X(k+N/2) +for k=1:(round(N/2)) + X1(k) = F1(k) + Wnk(k-1)*F2(k); + X2(k) = F1(k) - Wnk(k-1)*F2(k); +end +X_without_fft = [X1,X2] +%X_without_fft = 1×8 complex +% 36.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i +%⋯ +X_mag = abs(X_without_fft); +X_phase = angle(X_without_fft); +f = 1:N; +figure("Name","FFT") +subplot(3,1,1) +stem(f,x,"filled","LineWidth",3) +title("Input sequence") +ylabel("x(n)") +xlabel("N") +grid +subplot(3,1,2) +stem(f,X_mag,"filled","LineWidth",3) +title("Magnitude spectrum") +ylabel("Magnitude") +xlabel("Frequency") +grid +subplot(3,1,3) +stem(f,X_phase,"filled","LineWidth",3) +title("Phase spectrum") +ylabel("Phase") +xlabel("Frequency") +grid diff --git a/FIR.m b/FIR.m new file mode 100644 index 0000000..30aaf37 --- /dev/null +++ b/FIR.m @@ -0,0 +1,58 @@ +N=7; +wc=1; +alpha=(N-1)/2; +n=0:1:N-1; +%Desired impulse response +for i=0:1:(N-1) + if(i~=alpha) + hd(i+1)=(sin(wc*(i-3))/(pi*(i-3))); + else + hd(i+1)=wc/pi; + end +end +subplot(5,1,1); +stem(n,hd); +xlabel('n'); +ylabel('hd(n)'); +title('Desired impulse response'); +wHm = zeros(1,7); +%Hamming window +for j=0: 6 + wHm(j+1)=0.54-(0.46*cos(2*pi*j/(N-1))); +end +subplot(5,1,2); +stem(n,wHm); +xlabel('n'); +ylabel('whm(n)'); +title('Hamming window'); +%Impulse response +hn=hd.*wHm; +subplot(5,1,3); +stem(n,hn); +xlabel('n'); +ylabel('h(n)'); +title('Impulse response'); +hw=[]; +for w=0:(1/pi):pi + t=0; + temp=0; + const=hn(alpha+1); + for m=0:1:((N-3)/2) + temp=temp+(2.*hn(m+1).*cos(w.*(alpha-m))); + end + temp=temp+const; + hw=[hw, temp]; +end +w=0:1/pi:pi; +%Magnitude response +subplot(5,1,4); +stem(w,hw); +xlabel('w'); +ylabel('|H(w)|'); +title('Magnitude response'); +phase=-w*alpha; +subplot(5,1,5); +stem(w,phase); +xlabel('w'); +ylabel('Phase'); +title('Phase response'); \ No newline at end of file diff --git a/FIR1.m b/FIR1.m new file mode 100644 index 0000000..40c9124 --- /dev/null +++ b/FIR1.m @@ -0,0 +1,49 @@ +wc=1; +N=51; +alpha=(N-1)/2; +n=0:1:N-1; +%Rectangular window +for j=0:1:(N-1) + wR(j+1)=1; +end +subplot(5,1,1); +stem(n,wR); +xlabel('n'); +ylabel('wR(n)'); +title('Rectangular window'); +%Bartlett window +for j=0:1:(N-1) + wBr(j+1)=1-(2.*abs(j-((N-1)/2)))/(N-1); +end; +subplot(5,1,2); +stem(n,wBr); +xlabel('n'); +ylabel('wBr(n)'); +title('Bartlett window'); +%Hamming window +for j=0:1:(N-1) + wHm(j+1)=0.54-(0.46*cos(2*pi*j/(N-1))); +end; +subplot(5,1,3); +stem(n,wHm); +xlabel('n'); +ylabel('whm(n)'); +title('Hamming window'); +%Hanning window +for j=0:1:(N-1) + wHn(j+1)=0.5-(0.5*cos(2*pi*j/(N-1))); +end; +subplot(5,1,4); +stem(n,wHn); +xlabel('n'); +ylabel('whn(n)'); +title('Hanning window'); +%Blackmann window +for j=0:1:(N-1) + wBl(j+1)=0.42-(0.5*(cos(2*pi*j/(N-1))))+(0.08*(cos(4*pi*j/(N-1)))); +end; +subplot(5,1,5); +stem(n,wBl); +xlabel('n'); +ylabel('wbl(n)'); +title('Blackmann window'); \ No newline at end of file diff --git a/IDFT.m b/IDFT.m new file mode 100644 index 0000000..d7c4b5b --- /dev/null +++ b/IDFT.m @@ -0,0 +1,25 @@ +X = [1,2,3,4]; +N = length(X); +x = zeros(1,N); +j = sqrt(-1); +disp("IDFT of X(K)") +%IDFT of X(K) +for n=1:N + for k=1:N + x(n) = (x(n) + X(k)*(exp(j * n * k *(2 *pi / N)))) / N + end +end +ny = (1:N); +x_mag = abs(x); +x_phase = angle(x); +figure("Name","IDFT") +subplot(2,1,1) +stem(ny,x_mag,"LineWidth",3) +xlabel("Frequency") +ylabel("Mag") +title("Magnitude plot") +subplot(2,1,2) +stem(ny, x_phase,"LineWidth",3) +xlabel("Frequency") +ylabel("Phase") +title("Phase plot") \ No newline at end of file diff --git a/IFFT.m b/IFFT.m new file mode 100644 index 0000000..bd58914 --- /dev/null +++ b/IFFT.m @@ -0,0 +1,73 @@ +disp("Inverse Fast Fourier Transform - DITIFFT") +%Inverse Fast Fourier Transform - DITIFFT +syms Wnkm(k,m) Wnk(k) +X = [11 1+j 0 1-2j 0 1+2j 0 1-j]; +X = conj(X); +N = length(X); +f1 = (zeros(1,round(N/2))); +f2 = (zeros(1,round(N/2))); +y=1; +%Decimating into EVEN and ODD +for i=1:N + if rem(i,2)==0 + f2(y) = X(i); + y=y+1; + else + f1(y) = X(i); + + end +end +j=sqrt(-1); +Wnkm(k,m) = exp((-j*4*pi*k*m)/N); +Wnk(k) = exp((-j*2*pi*k)/N); +F1 = (zeros(1,round(N/2))); +F2 = (zeros(1,round(N/2))); +%F1(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F1(k) = F1(k) + f1(m)*Wnkm(k-1,m-1); + end +end +%F2(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F2(k) = F2(k) + f2(m)*Wnkm(k-1,m-1); + end +end +x1 = (zeros(1,round(N/2))); +x2 = (zeros(1,round(N/2))); +%X1(k) = X(k) and X2(k) = X(k+N/2) +for k=1:(round(N/2)) + x1(k) = (F1(k) + Wnk(k-1)*F2(k))/N; + x2(k) = (F1(k) - Wnk(k-1)*F2(k))/N; +end +x_without_ifft = [x1,x2] +%x_without_ifft = 1×8 +% 1.8750 1.5518 0.6250 1.5518 0.8750 1.1982 2.1250 ⋯ +X_mag = abs(X); +X_phase = angle(X); +f = 1:N; +%FFT with using inbuilt function +X = [11 1+j 0 1-2j 0 1+2j 0 1-j]; +x_with_ifft = ifft(X) +%x_with_ifft = 1×8 + %1.8750 1.5518 0.6250 1.5518 0.8750 1.1982 2.1250 ⋯ +figure("Name","IFFT") +subplot(3,1,1) +stem(f,X_mag,"filled","LineWidth",3) +title("X(k) Magnitude spectrum") +ylabel("Magnitude") +xlabel("Frequency") +grid +subplot(3,1,2) +stem(f,X_phase,"filled","LineWidth",3) +title("X(k) Phase spectrum") +ylabel("Phase") +xlabel("Frequency") +grid +subplot(3,1,3) +stem(f,x_without_ifft,"filled","LineWidth",3) +title("Output sequence") +ylabel("x(n)") +xlabel("N") +grid diff --git a/IIIR.m b/IIIR.m new file mode 100644 index 0000000..f400f26 --- /dev/null +++ b/IIIR.m @@ -0,0 +1,14 @@ +t=linspace(0,1,1000); +x=sin(2*pi*t); +z=awgn(x,1); +N=40; +f=10; +fs=200; +w=2*pi*f/fs; +[b,a]=butter(N,w); +fvtool(b,a); +iir=filter(b,a,z); +subplot(2,1,1) +plot(t,z); +subplot(2,1,2) +plot(t,iir); \ No newline at end of file diff --git a/IIR.m b/IIR.m new file mode 100644 index 0000000..9a5cd75 --- /dev/null +++ b/IIR.m @@ -0,0 +1,15 @@ +t=linspace(0,1,1000); +x=sin(2*pi*t); +z=awgn(x,1); +N=40; +f=10; +fs=200; +w=2*pi*f/fs; +[b,a]=butter(N,w); +fvtool(b,a); +iir=filter(b,a,z); +subplot(2,1,1) +plot(t,z); +subplot(2,1,2) +plot(t,iir); + diff --git a/IIR_basic.m b/IIR_basic.m new file mode 100644 index 0000000..63aa851 --- /dev/null +++ b/IIR_basic.m @@ -0,0 +1,15 @@ +kp=-1; +ks=-20; +wp=4; +ws=8; +n=(log((10^(-kp/10)-1)/(10^(-ks/10)-1)))/(2*(log(wp/ws))) +wc=(wp)/((10^(-kp/10)-1)^(1/10)) + +for k=(2*N)-1 + theta_k=((pi*k)/N)+(pi/(2*N))+(pi/2) +end +s0=((pi/(2*N))+(pi/2)) +s1=((pi*1)/N)+(pi/(2*N))+(pi/2) +s2=((pi*2)/N)+(pi/(2*N))+(pi/2) +s3=((pi*3)/N)+(pi/(2*N))+(pi/2) +s4=((pi*4)/N)+(pi/(2*N))+(pi/2) \ No newline at end of file diff --git a/SP_OPEN_ENDED.m b/SP_OPEN_ENDED.m new file mode 100644 index 0000000..07b2f9d --- /dev/null +++ b/SP_OPEN_ENDED.m @@ -0,0 +1,73 @@ +BW = "1 octave"; % Bandwidth +N = 8; % Filter order +F0 = 1000; % Center frequency (Hz) +Fs = 48000; % Sampling frequency (Hz) +of = octaveFilter(FilterOrder=N,CenterFrequency=F0, ... + Bandwidth=BW,SampleRate=Fs); +Nx = 100000; +scope1 = spectrumAnalyzer(SampleRate=Fs,Method="filter-bank", ... + AveragingMethod="exponential",PlotAsTwoSidedSpectrum=false, ... + FrequencyScale="log",FrequencySpan="start-and-stop-frequencies", ... + StartFrequency=1,StopFrequency=Fs/2,YLimits=[-60 10], ... + RBWSource="property",RBW=1); +tic +while toc < 20 + % Run for 20 seconds + x1 = randn(Nx,1); + y1 = of(x1); + scope1(y1) +end +ofb = octaveFilterBank("1/3 octave",Fs,FilterOrder=N); +freqz(ofb,NFFT=2^16) % Increase FFT length for better low-frequency resolution +set(gca,XScale="log") +axis([20 Fs/2 -50 5]) +title("1/3-Octave Filter Bank Magnitude Response") + +pinkNoise = dsp.ColoredNoise(Color="pink", ... + SamplesPerFrame=Nx, ... + NumChannels=1); + +scope2 = spectrumAnalyzer(SampleRate=Fs,Method="filter-bank", ... + AveragingMethod="exponential",PlotAsTwoSidedSpectrum=false, ... + FrequencyScale="log",FrequencySpan="start-and-stop-frequencies", ... + StartFrequency=20,StopFrequency=Fs/2,YLimits=[-40 30], ... + RBWSource="property",RBW=10); + +centerOct = getCenterFrequencies(ofb); +nbOct = numel(centerOct); +bandPower = zeros(1,nbOct); +nbSamples = 0; + +tic +while toc < 10 + xp = pinkNoise(); + yp = ofb(xp); + bandPower = bandPower + sum(yp.^2,1); + nbSamples = nbSamples + Nx; + scope2(yp) +end +b = 10^(3/10); % base-10 octave ratio +% Compute power (including pressure reference) +octPower = 10*log10(bandPower/nbSamples/4e-10); + +bar(log(centerOct)/log(b),octPower); +set(gca,Xticklabel=round(b.^get(gca,"Xtick"),2,"significant")); +title("1/3-Octave Power Spectrum") +xlabel("Octave Frequency Band (Hz)") +ylabel("Power (dB)") + +spl = splMeter(Bandwidth="1/3 octave", ... + OctaveFilterOrder=N, ... + SampleRate=Fs, ... + FrequencyWeighting="z-weighting"); + +scope3 = dsp.ArrayPlot(Title="Pink Noise SPL", ... + XLabel="Octave Frequency Band Number", ... + YLabel="Power (dB)",YLimits=[0 100]); +tic +while toc < 10 + xp = pinkNoise(); + yp = spl(xp); + ypm = mean(yp,1).'; + scope3(ypm) +end \ No newline at end of file diff --git a/SSPP.m b/SSPP.m new file mode 100644 index 0000000..3c1cc95 --- /dev/null +++ b/SSPP.m @@ -0,0 +1,72 @@ +N = 99; +[LPAnalysis, HPAnalysis, LPSynthsis, HPSynthesis] = firpr2chfb(N, 0.45); +fvt = fvtool(LPAnalysis,1, HPAnalysis,1, LPSynthsis,1, HPSynthesis,1); +fvt.Color = [1,1,1]; +legend(fvt,'Hlp Lowpass Decimator','Hhp Highpass Decimator',... + 'Glp Lowpass Interpolator','Ghp Highpass Interpolator'); +% Analysis section +analysisFilter = dsp.SubbandAnalysisFilter(LPAnalysis, HPAnalysis); +% Synthesis section +synthFilter = dsp.SubbandSynthesisFilter(LPSynthsis, HPSynthesis); +x = zeros(50,1); +x(1:3) = 1; +x(8:10) = 2; +x(16:18) = 3; +x(24:26) = 4; +x(32:34) = 3; +x(40:42) = 2; +x(48:50) = 1; +sigsource = dsp.SignalSource('SignalEndAction', 'Cyclic repetition',... + 'SamplesPerFrame', 50); +sigsource.Signal = x; + +% Scope to compare input signal with reconstructed output +sigcompare = dsp.ArrayPlot('NumInputPorts', 2, 'ShowLegend', true,... + 'Title', 'Input and reconstructed signals',... + 'ChannelNames',{'Input signal','Reconstructed signal'}); + +% Scope to plot the RMS error between the input and reconstructed signals +errorPlot = timescope('Title', 'RMS Error', 'SampleRate', 1, ... + 'TimeUnits', 'seconds', 'YLimits', [-0.5 2],... + 'TimeSpanSource', 'property','TimeSpan',100,... + 'TimeSpanOverrunAction','scroll'); + +% To calculate the transfer function of the cascade of Analysis and +% Synthesis subband filters +tfestimate = dsp.TransferFunctionEstimator('FrequencyRange','centered',... + 'SpectralAverages', 50); +% Scope to plot the magnitude response of the estimated transfer function +magplot = dsp.ArrayPlot('PlotType','Line', ... + 'YLabel', 'Magnitude Response (dB)',... + 'Title','Magnitude response of the estimated transfer function',... + 'XOffset',-25, 'XLabel','Frequency (Hz)',... + 'YLimits',[-5 5]); +% Scope to plot the phase response of the estimated transfer function +phaseplot = dsp.ArrayPlot('PlotType','Line', ... + 'YLabel', 'Phase Response (degrees)',... + 'Title','Phase response of the estimated transfer function',... + 'XOffset',-25, 'XLabel','Frequency (Hz)',... + 'YLimits',[-180 180]); + +for i = 1:100 + % Use the same signal as input in each iteration. + input = sigsource(); + % Analysis + [hi, lo] = analysisFilter(input); + % Synthesis + reconstructed = synthFilter(hi, lo); + + % Compensate for the delay caused by the filters and compare the + % signals. Since input signal is periodic, compare the current + % frames of input and output signals. + sigcompare(input(2:end), reconstructed(1:end-1)); + + % Plot error between signals + err = rms(input(2:end) - reconstructed(1:end-1)); + errorPlot(err); + + % Estimate transfer function of cascade + Txy = tfestimate(input(2:end), reconstructed(1:end-1)); + magplot(20*log10(abs(Txy))); + phaseplot(angle(Txy)*180/pi); +end diff --git a/circular_convolution.m b/circular_convolution.m new file mode 100644 index 0000000..d1682a0 --- /dev/null +++ b/circular_convolution.m @@ -0,0 +1,23 @@ +x = [1,2,3]; +h = [3,4,5,6]; +n1=length(x); +n2 = length(h); +n = max(n1,n2); +a = 1:n; +x = [x,zeros(1,n-n1)]; +h = [h,zeros(1,n-n2)]; +y = zeros(1,n); +for i =0:n-1 + for j = 0:n-1 + k = mod((i-j),n); + y(i+1) = y(i+1) + x(j+1)*h(k+1); + end +end +figure("Name","Circular conv"); +subplot(2,1,1); +stem(a,y,"LineWidth",3) +title("Circular Convolution: Without Inbuilt Function") +y = cconv(x,h,6); +subplot(2,1,2); +stem(y,"LineWidth",3) +title("Circular Convolution: With Inbuilt Function"); diff --git a/comb_func.m b/comb_func.m new file mode 100644 index 0000000..0d9224f --- /dev/null +++ b/comb_func.m @@ -0,0 +1,57 @@ +%dt comb +stem(ones(1,10), 'LineWidth', 2) +xlabel("Time"); +ylabel("Amplitude"); + +t = (-1: .1 : 1)'; +impulse = t == 0; +unitStep = t >= 0; +ramp = t .* unitStep; +exponential_inc = exp(1 * t); +exponential_dec = exp(-1 * t); +xlabel("Time"); +ylabel("Amplitude"); + +% Impulse unitStep Ramp +stem(t,[impulse unitStep ramp],"LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +hold off + +% Exponential +stem(t,exponential_inc,"LineWidth",3) +hold on +stem(t,exponential_dec,"LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +hold off + +% Rectangular pulse +fs = 1000; +T = 10*(1/50); +t = -0.1:.01:T-1/fs; +w = 20e-3; +x = rectpuls(t,w); +stem(t,x,"LineWidth",3) +xlabel("Time") +ylabel("Amplitude") +title("Rectangular") +grid on + +%Triangle +t = (0:0.01:0.4); +f = 10; +stem(t,sawtooth(2 * pi * f * t,1/2),"k-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); + +% Sinc And sa +x = linspace(-5,5,50); +sc = sinc(x); +sa = sin(x) ./ x; +stem(x,sc,"k-","LineWidth",3) +hold on +stem(x,sa,"r-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +hold off diff --git a/linear_circular.m b/linear_circular.m new file mode 100644 index 0000000..d9f954a --- /dev/null +++ b/linear_circular.m @@ -0,0 +1,14 @@ +x = [2 1 2 1]; +y = [1 2 3]; +clin = conv(x,y); +xpad = [x zeros(1,6-length(x))]; +ypad = [y zeros(1,6-length(y))]; +ccirc = ifft(fft(xpad).*fft(ypad)); +subplot(2,1,1) +stem(clin,'filled',"LineWidth",3) +ylim([0 11]) +title('Linear Convolution of x and y') +subplot(2,1,2) +stem(ccirc,'filled',"LineWidth",3) +ylim([0 11]) +title('Circular Convolution of xpad and ypad') \ No newline at end of file diff --git a/linear_conv_fft.m b/linear_conv_fft.m new file mode 100644 index 0000000..bd7f95c --- /dev/null +++ b/linear_conv_fft.m @@ -0,0 +1,53 @@ +x = [1 1 1 1]; +h = [1 0 1 0]; +n1 = length(x); +n2 = length(h); +N = n1+n2-1; +x = [x,zeros(1,N-n1)]; +h = [h,zeros(1,N-n2)]; +%To find X(k) and H(k) +X = fft(x); +H = fft(h); +%Y(k) +Y = X.*H; +%IFFT to get y(n) +syms Wnkm(k,m) Wnk(k) +Y = conj(Y); +f1 = (zeros(1,round(N/2))); +f2 = (zeros(1,round(N/2))); +y=1; +%Decimating into EVEN and ODD +for i=1:N + if rem(i,2)==0 + f2(y) = Y(i); + y=y+1; + else + f1(y) = Y(i); + + end +end +j=sqrt(-1); +Wnkm(k,m) = exp((-j*4*pi*k*m)/N); +Wnk(k) = exp((-j*2*pi*k)/N); +F1 = (zeros(1,round(N/2))); +F2 = (zeros(1,round(N/2))); +%F1(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F1(k) = F1(k) + f1(m)*Wnkm(k-1,m-1); + end +end +%F2(k) +for k=1:(round(N/2)) + for m=1:(round(N/2)) + F2(k) = F2(k) + f2(m)*Wnkm(k-1,m-1); + end +end +y1 = (zeros(1,round(N/2))); +y2 = (zeros(1,round(N/2))); +%X1(k) = X(k) and X2(k) = X(k+N/2) +for k=1:(round(N/2)) + y1(k) = (F1(k) + Wnk(k-1)*F2(k))/N; + y2(k) = (F1(k) - Wnk(k-1)*F2(k))/N; +end +y_linear_convolved = [y1,y2] \ No newline at end of file diff --git a/linear_convolution.m b/linear_convolution.m new file mode 100644 index 0000000..38a6e3f --- /dev/null +++ b/linear_convolution.m @@ -0,0 +1,41 @@ +x = [1,2,3]; +h = [3,4,5,6]; +x_length = length(x); +h_length = length(h); +N = x_length + h_length -1; +x = [x,zeros(1,N-x_length)]; +h = [h,zeros(1,N-h_length)]; +y = zeros(1,N); + +%Without inbuilt function +for n=1:N +for m=1:n +y(n) = y(n)+x(m)*h(n-m+1); +end +end +figure("Name","without-Inbuilt") +ny = (0:N-1); +subplot(2,2,1); +stem(ny,x,"LineWidth",3); +xlabel("n"); +ylabel("x(n)"); +title("First sequence"); +subplot(2,2,2); +stem(ny,h,"LineWidth",3); +xlabel("n"); +ylabel("h(n)"); +title("Second sequence"); +subplot(2,2,3); +stem(ny,y,"LineWidth",3); +xlabel("n"); +ylabel("y(n)"); +title("Convoluted signal"); + +%With inbuilt function +x = [1,2,3]; +h = [3,4,5,6]; +y = conv(x,h); +ny = (0:N-1); +subplot(2,2,4) +stem(ny,y,"LineWidth",3); +title('Inbuilt Function') diff --git a/sampling_func.m b/sampling_func.m new file mode 100644 index 0000000..b9c99cb --- /dev/null +++ b/sampling_func.m @@ -0,0 +1,49 @@ +%To Demonstrate sampling theorm +f=2; +fs1=1.2*f; +fs2=2*f; +fs3=8*f; +t=0:0.01:3; +t1=0:1/fs1:3; +t2=0:1/fs2:3; +t3=0:1/fs3:3; +x=sin(2*pi*f*t); +subplot(2,2,1); +plot(t,x); +xlabel('Time'); +ylabel('Amplitude'); +title('Continuous sine wave'); + +%Under sampling +% for fs<2fm +y1=sin(2*pi*f*t1); +subplot(2,2,2); +stem(t1,y1); +hold on +plot(t1,y1) +xlabel('Time'); +ylabel('Amplitude'); +title('Sampled output fs<2fm'); +hold off + +%Critical sampling for fs=2fm +y2=sin(2*pi*f*t2); +subplot(2,2,3); +stem(t2,y2); +hold on +plot(t2,y2); +xlabel('Time'); +ylabel('Amplitude'); +title('Sampled output fs=2fm'); +hold off + +%Over sampling for fs>2fm +y3=sin(2*pi*f*t3); +subplot(2,2,4); +stem(t3,y3); +hold on +plot(t3,y3); +xlabel('Time'); +ylabel('Amplitude'); +title('Sampled output fs>2fm'); +hold off \ No newline at end of file diff --git a/signals.m b/signals.m new file mode 100644 index 0000000..765e810 --- /dev/null +++ b/signals.m @@ -0,0 +1,54 @@ +t = (-1: .01 : 1)'; +impulse = t == 0; +unitStep = t >= 0; +ramp = t .* unitStep; +exponential_inc = exp(1 * t); +exponential_dec = exp(-1 * t); + +% Impulse unitStep Ramp +plot(t,[impulse unitStep ramp],"LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +legend("Impulse","Unit step","Ramp") +hold off + +% Exponential +plot(t,exponential_inc,"LineWidth",3) +hold on +plot(t,exponential_dec,"LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +legend("Rise","Fall") +hold off + +% Sinc And sa +x = linspace(-5,5); +sc = sinc(x); +sa = sin(x) ./ x; +plot(x,sc,"k-","LineWidth",3) +hold on +plot(x,sa,"r-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +legend("sinc","sa") +hold off +grid on + +% Rectangular pulse +fs = 1000; +T = 10*(1/50); +t = -0.1:.01:T-1/fs; +w = 20e-3; +x = rectpuls(t,w); +plot(t,x,"LineWidth",3) +xlabel("Time") +ylabel("Amplitude") +title("Rectangular") +grid on + +%Triangle +t = (0:0.001:0.4); +f = 10; +plot(t,sawtooth(2 * pi * f * t,1/2),"k-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); \ No newline at end of file diff --git a/sine_func.m b/sine_func.m new file mode 100644 index 0000000..59fb327 --- /dev/null +++ b/sine_func.m @@ -0,0 +1,15 @@ +t = (0:.01:2); +f = 1; +plot(t,sin(2 * pi * f * t),"k-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +grid on + +%ct_dt sine +fs = 40; +t = linspace(0,2,fs*2); +f = 1; +stem(t,sin(2 * pi * f * t),"k-","LineWidth",3) +xlabel("Time"); +ylabel("Amplitude"); +