Commot DSP

DSP/Synthetic_design/fft_Hypo.m

+ function  [k,xk]=fft_Hypo(fs,N,data) %FFT函数，已去直流分量
+ data=data-mean(data);
+ data_fft=fft(data,N);%fft变换
+ k=((1:N/2)-1)*fs/N;
+ xk=abs(data_fft(1:N/2));
\ No newline at end of file

DSP/Synthetic_design/filter_FIR_ButterBPF.m

+function Hd = filter_FIR_ButterBPF(Fs,Fstop1,Fpass1,Fpass2,Fstop2)
+%FILTER_FIR_BUTTERBPF Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 8.6 and the Signal Processing Toolbox 7.1.
+% Generated on: 03-Dec-2017 21:39:40
+
+% FIR Window Bandpass filter designed using the FIR1 function.
+
+% All frequency values are in Hz.
+%Fs = 1000;  % Sampling Frequency
+
+%Fstop1 = 50;              % First Stopband Frequency
+%Fpass1 = 90;              % First Passband Frequency
+%Fpass2 = 110;             % Second Passband Frequency
+%Fstop2 = 150;             % Second Stopband Frequency
+Dstop1 = 0.001;           % First Stopband Attenuation
+Dpass  = 0.057501127785;  % Passband Ripple
+Dstop2 = 0.0001;          % Second Stopband Attenuation
+flag   = 'scale';         % Sampling Flag
+
+% Calculate the order from the parameters using KAISERORD.
+[N,Wn,BETA,TYPE] = kaiserord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 ...
+                             1 0], [Dstop1 Dpass Dstop2]);
+
+% Calculate the coefficients using the FIR1 function.
+b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
+Hd = dfilt.dffir(b);
+
+% [EOF]

DSP/Synthetic_design/filter_IIR_Butter.m

+function Hd = filter_IIR_Butter(Fs,N,Fc1,Fc2)
+%FILTER_IIR_BUTTER Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 8.6 and the DSP System Toolbox 9.1.
+% Generated on: 30-Oct-2017 18:25:15
+
+% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.
+
+% All frequency values are in Hz.
+
+%%%%%%%Fs = 1000;  % Sampling Frequency
+
+%%%%%%%%N   = 8;   % Order
+%%%%%%%%Fc1 = 4;   % First Cutoff Frequency
+%%%%%%%%Fc2 = 15;  % Second Cutoff Frequency
+
+% Construct an FDESIGN object and call its BUTTER method.
+h  = fdesign.bandpass('N,F3dB1,F3dB2', N, Fc1, Fc2, Fs);
+Hd = design(h, 'butter');
+
+% [EOF]

DSP/Synthetic_design/filter_butterLPF.m

+function Hd = filter_butterLPF(Fs,N,Fc)
+%FILTER_LOWPASS Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 8.6 and the DSP System Toolbox 9.1.
+% Generated on: 09-Nov-2017 21:00:59
+
+% Butterworth Lowpass filter designed using FDESIGN.LOWPASS.
+
+% All frequency values are in Hz.
+%%%%%%%%%Fs = 10000;  % Sampling Frequency
+
+%%%%%%%%%%N  = 10;   % Order
+%%%%%%%%%%Fc = 100;  % Cutoff Frequency
+
+% Construct an FDESIGN object and call its BUTTER method.
+h  = fdesign.lowpass('N,F3dB', N, Fc, Fs);
+Hd = design(h, 'butter');
+
+% [EOF]

DSP/Synthetic_design/main.m

+ function  homework(fs, N) %fs=1000;N=200       homework(1000,400) 
+ 
+n=0:N-1;
+t=n/fs;
+figure(1);
+base1=0.5*square(2*pi*10*t)+0.3;subplot(3,4,1);plot(t,base1);title('基波1(10HZ方波)')
+carrier1=sin(2*pi*100*t);subplot(3,4,2);plot(t,carrier1);title('载波1(100HZ正弦波)')
+mix1=base1.*carrier1;subplot(3,4,3);plot(t,mix1);title('基波1调制结果')
+
+%base2=sawtooth(2*pi*10*t,0.5)+1;subplot(3,4,5);plot(t,base2);title('基波10HZ三角波')
+base2=0.5*square(2*pi*10*t)+0.3;subplot(3,4,5);plot(t,base2);title('基波2（10HZ方波）')
+carrier2=sin(2*pi*300*t);subplot(3,4,6);plot(t,carrier2);title('载波2（300HZ正弦波）')
+mix2=base2.*carrier2;subplot(3,4,7);plot(t,mix2);title('基波2调制结果')
+
+mix3=mix1+mix2;subplot(3,4,8);plot(t,mix3);title('调制信号叠加结果')
+
+bandpass1=filter(filter_IIR_Butter(fs,4,50,150),mix3);subplot(3,4,9);plot(t,bandpass1);title('50~150HZ带通滤波');%IIR方式带通滤波
+%bandpass1=filter(filter_FIR_ButterBPF(fs,10,50,150,200),mix3);subplot(3,4,9);plot(t,bandpass1);title('50~150HZ带通滤波');%FIR方式带通滤波
+demod1=bandpass1.*carrier1;%相干解调
+final1=filter(filter_butterLPF(fs,12,100),demod1);subplot(3,4,10);plot(t,final1);title('100HZ低通&解调信号1“最终结果”');
+
+bandpass2=filter(filter_IIR_Butter(fs,4,250,350),mix3);subplot(3,4,11);plot(t,bandpass2);title('250~350HZ带通滤波');%IIR方式带通滤波
+%bandpass2=filter(filter_FIR_ButterBPF(fs,150,220,380,450),mix3);subplot(3,4,11);plot(t,bandpass2);title('220~380HZ带通滤波');%FIR方式带通滤波
+demod2=bandpass2.*carrier2;%相干解调
+final2=filter(filter_butterLPF(fs,12,100),demod2);subplot(3,4,12);plot(t,final2);title('100HZ低通&解调信号2“最终结果”');
+
+figure(2);
+[k,xk]=fft_Hypo(fs,N,base1);subplot(2,2,1);stem(k,xk);title('基波1频谱图(10HZ方波)');
+[k,xk]=fft_Hypo(fs,N,final1);subplot(2,2,2);stem(k,xk);title('解调信号1频谱图(10HZ方波)');
+[k,xk]=fft_Hypo(fs,N,base2);subplot(2,2,3);stem(k,xk);title('基波2频谱图(10HZ方波)');
+[k,xk]=fft_Hypo(fs,N,final2);subplot(2,2,4);stem(k,xk);title('基波2频谱图(10HZ方波)');
+figure(3);
+[k,xk]=fft_Hypo(fs,N,mix1);subplot(2,3,1);stem(k,xk);title('基波1调制结果频谱图');
+[k,xk]=fft_Hypo(fs,N,mix2);subplot(2,3,2);stem(k,xk);title('基波2调制结果频谱图');
+[k,xk]=fft_Hypo(fs,N,mix3);subplot(2,3,3);stem(k,xk);title('叠加结果频谱图');
+[k,xk]=fft_Hypo(fs,N,bandpass1);subplot(2,3,4);stem(k,xk);title('50~150HZ带通滤波频谱图');
+[k,xk]=fft_Hypo(fs,N,bandpass2);subplot(2,3,5);stem(k,xk);title('220~380HZ带通滤波频谱图');
+[k,xk]=fft_Hypo(fs,N,base1);subplot(2,3,6);stem(k,xk);title('测试');
+figure(4);
+[k,xk]=fft_Hypo(fs,N,demod1);subplot(1,2,1);stem(k,xk);title('信号1，相干解调后，低通前');
+[k,xk]=fft_Hypo(fs,N,demod2);subplot(1,2,2);stem(k,xk);title('信号2，相干解调后，低通前');
+%figure(5);
+%[k,xk]=fft_Hypo(fs,N,carrier1);stem(k,xk);title('测试');
+%fvtool(filter_FIR_ButterBPF(fs,50,80,120,150));
+%fvtool(filter_FIR_ButterBPF(fs,250,280,320,350));
+%fvtool(filter_butterLPF(fs,10,100));

DSP/experiment/impseq.m

+ function [x,n] = impseq(n0,n1,n2)
+ %产生x(n)=ε(n-n0) 时域范围n1 n2
+    % 产生 x(n) = delta(n-n0); n1 <= n,n0 <= n2
+    % [x,n] = impseq(n0,n1,n2)
+    if ((n0 < n1) | (n0 > n2) | (n1 > n2))
+        error('参数必须满足 n1 <= n0 <= n4')
+    end
+    n = [n1:n2];
+    %x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))];
+x = [(n-n0) == 0]; 

DSP/experiment/sigFourierTran.m

+ function [X,magX,angX] = sigFourierTran(x,n,dot)
+    % 计算离散序列的付立叶变换
+    % [X,magX,angX] = FourierTran(x,n)
+    % 或[X,magX,angX] = FourierTran(x,n,dot)
+    if nargin< 3   % nargin是输入变量的个数
+        dot=600;
+    end
+    k=-dot:dot;
+    w=(pi/dot)*k; 
+    X=x*(exp(-j).^(n'*w));%做傅里叶变换
+    magX=abs(X);
+    angX=angle(X);
+    subplot(211);
+    plot(w/pi,magX);
+    xlabel('频率(单位{\pi})');ylabel('|X(e^{ j\omega})|');
+    title('幅频特性');
+    subplot(212);
+    plot(w/pi,angX/pi);
+    xlabel('频率(单位{\pi})');ylabel('弧度/{\pi}');
+    title('相频特性'); 

DSP/experiment/sigadd.m

+ function [y,n] = sigadd(x1,n1,x2,n2)
+    % 实现 y(n) = x1(n)+x2(n)
+    % [y,n] = sigadd(x1,n1,x2,n2)
+    % y = 在包含n1 和 n2 的n点上求序列和 
+    % x1 = 在 n1上的第一序列
+% x2 = 在 n2上的第二序列(n2可与 n1不等)
+n = min(min(n1),min(n2)):max(max(n1),max(n2));  % y(n)的长度                                              
+y1 = zeros(1,length(n)); y2 = y1;               % 初始化
+y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;  % 具有y的长度的 x1                                                                  
+ y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;  % 具有y的长度的x2                                                                 
+y = y1+y2;                         % 序列相加.

DSP/experiment/sigconjsym.m

+function [y,n] = sigconjsym(x,n,type)%type=0求共轭对称分量，type=1求共轭反对称分量
+if(type==0)   %Xf=0.5*(X(n) + X*(-n))
+    [x1,n1]=sigfold(conj(x),n);
+    [y,n]=sigadd(x,n,x1,n1);
+    y=y/2;
+end 
+if(type==1)   %Xf=0.5*(X(n) - X*(-n))
+    [x1,n1]=sigfold(conj(x),n);
+    [y,n]=sigadd(x,n,-x1,n1);
+    y=y/2;
+end 
\ No newline at end of file

DSP/experiment/sigfold.m

+function [y,n] = sigfold(x,n)
+    % ʵ y(n) = x(-n)
+    % [y,n] = sigfold(x,n)
+    y = fliplr(x); n = -fliplr(n); 

DSP/experiment/sigmult.m

+function [y,n] = sigmult(x1,n1,x2,n2)
+    % 实现 y(n) = x1(n)*x2(n)
+    % [y,n] = sigmult(x1,n1,x2,n2)
+    % y = 在n区间上的乘积序列,n 包含 n1 和 n2
+    % x1 = 在 n1上的第一序列
+    % x2 = 在 n2上的第二序列(n2可与 n1不等)
+n = min(min(n1),min(n2)):max(max(n1),max(n2)); % y(n)的长度
+    y1 = zeros(1,length(n)); y2 = y1;                %初始化
+    y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;     % 具有y的长度的 x1 
+    y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;     % 具有y的长度的 x2 
+    y = y1 .* y2;                                  % 序列相乘 

DSP/experiment/sigscale.m

+function [y,n] = sigscale(x,n,m)%y(n)=x(mn)  横坐标缩小n倍  使用条件：m>0
+n1=n/m;
+n2=find(n1==fix(n1));
+n3=n1(n2);%找到变换后的n
+x3=x(n2);%找到变换后的x
+y=x3;n=n3;
\ No newline at end of file

DSP/experiment/sigshift.m

+function [y,n] = sigshift(x,m,n0) %将x(n0)向右平移m个单位得到y(n)
+    % 实现 y(n) = x(n-n0)
+    % [y,n] = sigshift(x,m,n0)
+n = m+n0; y = x;

DSP/experiment/stepseq.m

+function [x,n] = stepseq(n0,n1,n2)
+    %  x(n) = u(n-n0);
+    %n1 <= n0 <= n2
+    % [x,n] = stepseq(n0,n1,n2)
+    if ((n0 < n1) | (n0 > n2) | (n1 > n2))
+        error(' n1 <= n0 <= n2')
+    end
+    n = [n1:n2];
+    %x = [zeros(1,(n0-n1)), ones(1,(n2-n0+1))];
+x = [(n-n0) >= 0];