irrdf.m

通信系统中基于MATLAB的关于非线性滤波器的仿真程序

function varargout = irrdf(varargin)
% IRRDF M-file for irrdf.fig
%      IRRDF, by itself, creates a new IRRDF or raises the existing
%      singleton*.
%
%      H = IRRDF returns the handle to a new IRRDF or the handle to
%      the existing singleton*.
%
%      IRRDF('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IRRDF.M with the given input arguments.
%
%      IRRDF('Property','Value',...) creates a new IRRDF or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before irrdf_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to irrdf_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help irrdf

% Last Modified by GUIDE v2.5 03-Jan-2008 01:29:32

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @irrdf_OpeningFcn, ...
                   'gui_OutputFcn',  @irrdf_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before irrdf is made visible.
function irrdf_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to irrdf (see VARARGIN)

% Choose default command line output for irrdf
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes irrdf wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = irrdf_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
figure(1);
[y,Fs,bits]=wavread('wang1.wav');
subplot(2,2,1);
plot(y);title('1 原始信号时域波形');grid
fp=600;fs=1600;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Mp=tan(wp/2);Ms=tan(ws/2);
Ap=1;As=100;
[N,Wn]=buttord(Mp,Ms,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[numz,denz]=bilinear(nums,dens,0.5);
ylp=filter(numz,denz,y);
subplot(2,2,2);
plot(ylp);title('2 低通输出');grid
fs=4800;fp=6000;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Mp=tan(wp/2);Ms=tan(ws/2);
Wp=1;Ws=Mp/Ms;
Ap=1;As=100;
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2hp(nums,dens,Ws);
[numz,denz]=bilinear(bt,at,0.5);
yhp=filter(numz,denz,y);
subplot(2,2,3);
plot(yhp);title('3 高通输出');grid
fp1=1200;fs1=600;fp2=3000;fs2=3600;
wp1=2*pi*fp1/Fs;ws1=2*pi*fs1/Fs;wp2=2*pi*fp2/Fs;ws2=2*pi*fs2/Fs;
Ap=1;As=32;
Mp1=tan(wp1/2);Ms1=tan(ws1/2);Mp2=tan(wp2/2);Ms2=tan(ws2/2);
Bw=Mp2-Mp1;P=Mp1*Mp2; S=Ms1*Ms2;
if P>S Ms=(Mp1*Mp2)/Ms2;else Ms=(Mp1*Mp2)/Ms1;end
Wp=1;Ws=(Mp1*Mp2-Ms.^2)/(Ms*Bw);
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2bp(nums,dens,sqrt(Mp1*Mp2),Bw);
[numz,denz]=bilinear(bt,at,0.5);
ybp=filter(numz,denz,y);
subplot(2,2,4);
plot(ybp);title('4 带通输出');grid

figure(2);
[y,Fs,bits]=wavread('wang1.wav');
YDFT=fft(y);  
subplot(2,2,1);
plot(abs(YDFT));title('1 原始信号幅度谱');grid
fp=600;fs=1600;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Mp=tan(wp/2);Ms=tan(ws/2);
Ap=1;As=100;
[N,Wn]=buttord(Mp,Ms,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[numz,denz]=bilinear(nums,dens,0.5);
ylp=filter(numz,denz,y);
YlpDFT=fft(ylp);
subplot(2,2,2);
plot(abs(YlpDFT));title('2 低通输出频率幅度响应');grid
fs=4800;fp=6000;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Mp=tan(wp/2);Ms=tan(ws/2);
Wp=1;Ws=Mp/Ms;
Ap=1;As=100;
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2hp(nums,dens,Ws);
[numz,denz]=bilinear(bt,at,0.5);
yhp=filter(numz,denz,y);
YhpDFT=fft(yhp);
subplot(2,2,3);
plot(abs(YhpDFT));title('3 高通输出频率幅度响应');grid
fp1=1200;fs1=600;fp2=3000;fs2=3600;
wp1=2*pi*fp1/Fs;ws1=2*pi*fs1/Fs;wp2=2*pi*fp2/Fs;ws2=2*pi*fs2/Fs;
Ap=1;As=32;
Mp1=tan(wp1/2);Ms1=tan(ws1/2);Mp2=tan(wp2/2);Ms2=tan(ws2/2);
Bw=Mp2-Mp1;P=Mp1*Mp2; S=Ms1*Ms2;
if P>S Ms=(Mp1*Mp2)/Ms2;else Ms=(Mp1*Mp2)/Ms1;end
Wp=1;Ws=(Mp1*Mp2-Ms.^2)/(Ms*Bw);
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2bp(nums,dens,sqrt(Mp1*Mp2),Bw);
[numz,denz]=bilinear(bt,at,0.5);
ybp=filter(numz,denz,y);
YbpDFT=fft(ybp);
subplot(2,2,4);
plot(abs(YbpDFT));title('4带通输出频率幅度响应');grid

figure(3);
fp=600;fs=1600;
wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;
Mp=tan(wp/2);Ms=tan(ws/2);
Ap=1;As=100;
[N,Wn]=buttord(Mp,Ms,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[numz,denz]=bilinear(nums,dens,0.5);
[Hlp,w]=freqz(numz,denz);
subplot(1,3,1);plot(w/pi,20*log10(abs(Hlp)));grid
title('1 低通滤波器的幅度响应');xlabel('归一化频率(w/pi)');ylabel('幅度谱');
fs=4800;fp=6000;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Mp=tan(wp/2);Ms=tan(ws/2);
Wp=1;Ws=Mp/Ms;
Ap=1;As=100;
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2hp(nums,dens,Ws);
[numz,denz]=bilinear(bt,at,0.5);
[Hhp,w]=freqz(numz,denz);
subplot(1,3,2);plot(w/pi,20*log10(abs(Hhp)));grid
title('2 高通滤波器的幅度响应');xlabel('归一化频率(w/pi)');ylabel('幅度谱');
fp1=1200;fs1=600;fp2=3000;fs2=3600;
wp1=2*pi*fp1/Fs;ws1=2*pi*fs1/Fs;wp2=2*pi*fp2/Fs;ws2=2*pi*fs2/Fs;
Ap=1;As=32;
Mp1=tan(wp1/2);Ms1=tan(ws1/2);Mp2=tan(wp2/2);Ms2=tan(ws2/2);
Bw=Mp2-Mp1;P=Mp1*Mp2; S=Ms1*Ms2;
if P>S Ms=(Mp1*Mp2)/Ms2; else Ms=(Mp1*Mp2)/Ms1;end
Wp=1;Ws=(Mp1*Mp2-Ms.^2)/(Ms*Bw);
[N,Wn]=buttord(Wp,Ws,Ap,As,'s');
[nums,dens]=butter(N,Wn,'s');
[bt,at]=lp2bp(nums,dens,sqrt(Mp1*Mp2),Bw);
[numz,denz]=bilinear(bt,at,0.5);
[Hbp,w]=freqz(numz,denz);
subplot(1,3,3);plot(w/pi,20*log10(abs(Hbp)));grid
title('3 带通滤波器的幅度响应');xlabel('归一化频率(w/pi)');ylabel('幅度谱');