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"); +