www.gusucode.com > 声音的处理有：LPC，FFT，共振峰，频谱源码程序 > siganlandsystemusingMatlab/SSUM/fseries/fseriesexpo.m
function varargout = fseriesexpo(varargin)
% This code is by Jim Squire <SquireJC@vmi.edu>
% It needs to be recoded to SSUM format.
%
% FSERIESEXPO M-file for fseriesexpo.fig
% FSERIESEXPO, by itself, creates a new FSERIESEXPO or raises the existing
% singleton*.
%
% H = FSERIESEXPO returns the handle to a new FSERIESEXPO or the handle to
% the existing singleton*.
%
% FSERIESEXPO('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in FSERIESEXPO.M with the given input arguments.
%
% FSERIESEXPO('Property','Value',...) creates a new FSERIESEXPO or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before fseriesexpo_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to fseriesexpo_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 fseriesexpo
% Last Modified by GUIDE v2.5 31-Mar-2005 18:04:20
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @fseriesexpo_OpeningFcn, ...
'gui_OutputFcn', @fseriesexpo_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin & isstr(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 fseriesexpo is made visible.
function fseriesexpo_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);
t=load('fseriespict.mat');
% Fig1
axes(handles.uiaxFig1);
image(t.strGraphic.Fig1);
set(handles.uiaxFig1,'Visible','Off')
colormap([0 0 0;236 233 216]/255)
% Fig2
axes(handles.uiaxFig2)
image(t.strGraphic.Fig2)
set(handles.uiaxFig2,'Visible','Off')
% Fig3
axes(handles.uiaxFig3)
image(t.strGraphic.Fig3)
set(handles.uiaxFig3,'Visible','Off')
% Calculate and fill the plots
CalculateAndUpdate(handles)
function varargout = fseriesexpo_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
%%%%%%%%%%%%%%%%%%%%%%%% Input function: wave type %%%%%%%%%%%%%%%%%%%%%%%%%%%
function uirbWaveTypeSin_Callback(hObject, eventdata, handles)
set(handles.uirbWaveTypeSin, 'Value', 1);
set([handles.uirbWaveTypeSquare handles.uirbWaveTypeTriangle], 'Value', 0);
CalculateAndUpdate(handles)
function uirbWaveTypeSquare_Callback(hObject, eventdata, handles)
set(handles.uirbWaveTypeSquare, 'Value', 1);
set([handles.uirbWaveTypeSin handles.uirbWaveTypeTriangle], 'Value', 0);
CalculateAndUpdate(handles)
function uirbWaveTypeTriangle_Callback(hObject, eventdata, handles)
set(handles.uirbWaveTypeTriangle, 'Value', 1);
set([handles.uirbWaveTypeSin handles.uirbWaveTypeSquare], 'Value', 0);
CalculateAndUpdate(handles)
%%%%%%%%%%%%%%%%%%%%%%% Slider data %%%%%%%%%%%%%%%%%%%%%%%%%
function uislFrequency_Callback(hObject, eventdata, handles)
CalculateAndUpdate(handles)
function uislAmplitude_Callback(hObject, eventdata, handles)
CalculateAndUpdate(handles)
function uislSymmetry_Callback(hObject, eventdata, handles)
CalculateAndUpdate(handles)
function uislYOffset_Callback(hObject, eventdata, handles)
CalculateAndUpdate(handles)
function uislTDelay_Callback(hObject, eventdata, handles)
CalculateAndUpdate(handles)
%%%%%%%%%%%%%%%%%%%%%%% Fourier series form %%%%%%%%%%%%%%%%%%%%%%%%%
function uirbFig1_Callback(hObject, eventdata, handles)
set(handles.uirbFig1, 'Value', 1)
set([handles.uirbFig2 handles.uirbFig3], 'Value', 0);
set(handles.uirbPlotComponent1,'String','ai')
set(handles.uirbPlotComponent2,'String','bi')
set([handles.uirbPlotComponent3 handles.uirbPlotComponent4],'Enable','off')
if get(handles.uirbPlotComponent3,'Value') ||get(handles.uirbPlotComponent4,'Value')
set([handles.uirbPlotComponent3 handles.uirbPlotComponent4],'Value',0)
set(handles.uirbPlotComponent1,'Value',1)
end
CalculateAndUpdate(handles)
function uirbFig2_Callback(hObject, eventdata, handles)
set(handles.uirbFig2, 'Value', 1);
set([handles.uirbFig1 handles.uirbFig3], 'Value', 0);
set(handles.uirbPlotComponent1,'String','Ai')
set(handles.uirbPlotComponent2,'String','Thetai')
set([handles.uirbPlotComponent3 handles.uirbPlotComponent4],'Enable','off')
if get(handles.uirbPlotComponent3,'Value') ||get(handles.uirbPlotComponent4,'Value')
set([handles.uirbPlotComponent3 handles.uirbPlotComponent4],'Value',0)
set(handles.uirbPlotComponent1,'Value',1)
end
CalculateAndUpdate(handles)
function uirbFig3_Callback(hObject, eventdata, handles)
set(handles.uirbFig3, 'Value', 1);
set([handles.uirbFig1 handles.uirbFig2], 'Value', 0);
set(handles.uirbPlotComponent1,'String','real{ci}')
set(handles.uirbPlotComponent2,'String','imag{ci}')
set([handles.uirbPlotComponent3 handles.uirbPlotComponent4],'Enable','on')
CalculateAndUpdate(handles)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fourier Series: NumTerms %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function uislNumTerms_Callback(hObject, eventdata, handles)
switch round(get(handles.uislNumTerms,'Value'))
case 0, N = 0;
case 1, N = 1;
case 2, N = 2;
case 3, N = 3;
case 4, N = 4;
case 5, N = 5;
case 6, N = 6;
case 7, N = 7;
case 8, N = 8;
case 9, N = 9;
case 10, N = 10;
case 11, N = 20;
case 12, N = 30;
case 13, N = 40;
case 14, N = 50;
case 15, N = 100;
end
set(handles.uitxNumTerms,'String',sprintf('N terms = %g',N))
CalculateAndUpdate(handles)
%%%%%%%%%%%%%%%%%%%%%%%%%% Fourier Series: plot components %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function uirbPlotComponent1_Callback(hObject, eventdata, handles)
set(handles.uirbPlotComponent1, 'Value', 1);
set([handles.uirbPlotComponent2 handles.uirbPlotComponent3 handles.uirbPlotComponent4], 'Value', 0);
CalculateAndUpdate(handles)
function uirbPlotComponent2_Callback(hObject, eventdata, handles)
set(handles.uirbPlotComponent2, 'Value', 1);
set([handles.uirbPlotComponent1 handles.uirbPlotComponent3 handles.uirbPlotComponent4], 'Value', 0);
CalculateAndUpdate(handles)
function uirbPlotComponent3_Callback(hObject, eventdata, handles)
set(handles.uirbPlotComponent3, 'Value', 1);
set([handles.uirbPlotComponent1 handles.uirbPlotComponent2 handles.uirbPlotComponent4], 'Value', 0);
CalculateAndUpdate(handles)
function uirbPlotComponent4_Callback(hObject, eventdata, handles)
set(handles.uirbPlotComponent4, 'Value', 1);
set([handles.uirbPlotComponent1 handles.uirbPlotComponent2 handles.uirbPlotComponent3], 'Value', 0);
CalculateAndUpdate(handles)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CalculateAndUpdate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CalculateAndUpdate(handles)
strGui = GetGUIValues(handles); % get all the scaled date from the GUI
% calculate strData.y0Time, .xTime, y1Time, xFreq, yFreq
strData = CaculateSeries(strGui);
% Display time
axes(handles.uiaxTime)
plot(strData.xTime,strData.y0Time,'r-',strData.xTime,strData.y1Time,'b-')
xlabel('Time (s)')
legend('Orignal','Reconstructed')
v=axis;
if v(4)-v(3)<1e-4
v(4)=v(4)+0.5;
v(3)=v(3)-0.5;
end
axis(v)
% Display freq
axes(handles.uiaxFreq)
stem('v6',strData.xFreq,strData.yFreq)
xlabel('Freq (Hz)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GetGUIValues %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function strGui = GetGUIValues(handles)
% Wave type
t = get(handles.uirbWaveTypeSin,'Value') + 2*get(handles.uirbWaveTypeSquare,'Value') + 4*get(handles.uirbWaveTypeTriangle,'Value');
switch t
case 1, strGui.waveType = 'Sin';
case 2, strGui.waveType = 'Square';
case 4, strGui.waveType = 'Triangle';
otherwise, error(sprintf('Unknown wave type: ',t));
end
% Wave parameters
strGui.w = ScaleIt(get(handles.uislFrequency,'Value'),0,4,8,0,2*pi,6*pi); % frequency in rads/s
if ~strGui.w, strGui.w=1e-6; end
T = 2*pi/strGui.w;
strGui.A = ScaleIt(get(handles.uislAmplitude,'Value'),0,4,8,0,1,4); % amplitude
strGui.sym = ScaleIt(round(get(handles.uislSymmetry,'Value')),0,4,8,1e-6,0.5-1e-6,1-1e-6); % symmetry
strGui.yOff = ScaleIt(get(handles.uislYOffset,'Value'),0,4,8,-2,0,2); % yOffset
strGui.tDelay = ScaleIt(get(handles.uislTDelay,'Value'),0,4,8,0,T/2,T); % tDelay
% Fourier series form
t = get(handles.uirbFig1,'Value') + 2*get(handles.uirbFig2,'Value') + 4*get(handles.uirbFig3,'Value');
switch t
case 1, strGui.fourierForm = 'Quadrature';
case 2, strGui.fourierForm = 'Polar';
case 4, strGui.fourierForm = 'Exponential';
otherwise, error(sprintf('Unknown Fourier Form',t));
end
% Num Terms
sNTerms = get(handles.uitxNumTerms,'String');
sNTerms(1:10) = [];
strGui.N = str2num(sNTerms);
% Plot component ai, bi, etc.
t = get(handles.uirbPlotComponent1,'Value') + 2*get(handles.uirbPlotComponent2,'Value') + 4*get(handles.uirbPlotComponent3,'Value') + 8*get(handles.uirbPlotComponent4,'Value');
switch t
case 1, strGui.plotComponent = 1;
case 2, strGui.plotComponent = 2;
case 4, strGui.plotComponent = 3;
case 8, strGui.plotComponent = 4;
otherwise, error(sprintf('Unknown component type %g',t))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CaculateSeries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function strData = CaculateSeries(strGui)
% get GUI parameters
waveType = strGui.waveType;
w = strGui.w;
T = 2*pi/w;
A = strGui.A;
sym = strGui.sym;
yOff = strGui.yOff;
tDelay = strGui.tDelay;
fourierForm = strGui.fourierForm;
N = strGui.N;
plotComponent = strGui.plotComponent;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Calculate the time-domain vectors %%%%%%%%%%%%%%%%%%%%%
% create output vectors
t = linspace(0,2,400);
y0 = zeros(size(t));
y1 = y0;
% First create a triangle wave, regardless of what was requested
T1 = T*sym; % Triangle changes direction at T1, which may be from 0 to T
m1 = -2/T1; b1 = 1; % line equations for first part
m2 = 2/(T-T1); b2 = 1-m2*T;
tt = mod(t-tDelay,T) + tDelay;
for i=1:length(t)
if tDelay<=tt(i) && tt(i)<=T1+tDelay
triangle(i) = m1*(tt(i)-tDelay) + b1;
elseif T1+tDelay < tt(i) && tt(i) <= T+tDelay
triangle(i) = m2*(tt(i)-tDelay) + b2;
else
error('Impossible time bracket')
end
end
% Next change it into the requested type
switch strGui.waveType
case 'Sin'
for i=1:length(y0)
theta = (1-triangle)*pi/2;
y0 = A*cos(theta) +yOff;
end
case 'Square'
if sym<=1e-6,
y0 = ones(size(triangle));
elseif sym>=1-1e-6
y0 = -ones(size(triangle));
else
symOffset = (sym-1/2)*2;
y0 = A*sign(triangle-symOffset) + yOff;
end
Ta = (m1*tDelay-b1+2*sym-1)/m1;
Tb = (m2*tDelay-b2+2*sym-1)/m2;
if w<=1e-6
Ta=0; Tb=0;
end
case 'Triangle'
y0 = A*triangle + yOff;
otherwise
error('Unknown waveType')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Calculate the c's %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create vectors
f = zeros(1,N+1);
c = f;
% determine c
switch strGui.waveType
case 'Sin'
tt1=m1*m2*pi*T*yOff;
tt2=m2*cos(b1*pi/2)-m1*cos(pi/2*(b2+m2*T));
tt3=-m2*cos(pi/2*(b1+m1*T1))+m1*cos(pi/2*(b2+m2*T1));
tt4=m1*m2*pi*T;
c(1)=(tt1+2*A*(tt2+tt3))/tt4;
if w<=1e-5
c(1)=y0(1);
else
for n=1:N
tt1=-1/(n*T*w);
tt2=pi/2;
tt3=m1^2*pi^2-4*n^2*w^2;
tt4=m2^2*pi^2-4*n^2*w^2;
tt5=tt2*(b1+m1*T1);
tt6=tt2*(b2+m2*T);
tt7=tt2*(b2+m2*T1);
tt8=j*m1*pi*cos(b1*tt2)-2*n*w*sin(b1*tt2);
tt9=-j*m1*pi*cos(tt5)+exp(j*n*T1*w)*tt8+2*n*w*sin(tt5);
tt10=-j*m2*pi*cos(tt6)+2*n*w*sin(tt6);
tt11=j*m2*pi*tt3*cos(tt7)+tt4*tt9-2*n*w*tt3*sin(tt7);
tt12=exp(-j*n*w*(T+tDelay))*yOff*(exp(j*n*T*w)-1);
tt13=exp(j*n*w*(T1+2*tDelay))*tt3*tt10+exp(j*n*w*(T+2*tDelay))*tt11;
tt14=2*A*n*w*tt13*exp(-j*n*w*(T+T1+3*tDelay));
c(n+1)=tt1*j*(tt14/(tt3*tt4));
end
end
case 'Triangle'
if w<=1e-5
c(1)=y0(1);
else
c(1)=A*(m2*T^2+2*b2*(T-T1)+2*b1*T1+(m1-m2)*T1^2)/(2*T)+yOff;
for n=1:N
tt1=j*n*w;
tt2=exp(-tt1*(T+T1+tDelay))/(n^2*T*w^2);
tt3=-exp(tt1*(T+T1))*(m1+tt1*b1);
tt4=exp(tt1*T1)*(m2+tt1*(b2+m2*T));
tt5=exp(tt1*T)*(m1-m2+tt1*(b1-b2+m1*T1-m2*T1));
tt6=-tt1*exp(tt1*T1)*(exp(tt1*T)-1)*yOff;
c(n+1)=tt2*(A*(tt3+tt4+tt5)+tt6);
end
end
case 'Square'
c(1)=A*(T+2*(Ta-Tb))/T+yOff; % this is really c0
for n=1:N
tt1=j*n*w;
tt2=j*exp(-tt1*(T+2*Ta+Tb))/(n*T*w);
tt3=-2*A*exp(tt1*(T+2*Ta));
tt4=A*exp(tt1*(Ta+Tb));
tt5=(A-yOff)*exp(tt1*(T+Ta+Tb));
tt6=yOff*exp(tt1*(Ta+Tb));
c(n+1)=tt2*(tt3+tt4+tt5+tt6);
end
otherwise
error('unknown waveType')
end
% if elements approximately zero, make exactly zero
for n=1:N+1
if abs(real(c(n)))<=1e-5
c(n)=j*imag(c(n));
end
if abs(imag(c(n)))<=1e-5
c(n)=real(c(n));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Synthesize f(t) %%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1=y1+c(1);
for n=1:N
y1=y1 + c(n+1)*exp(j*n*w*t) + conj(c(n+1))*exp(-j*n*w*t);
end
switch fourierForm
case 'Quadrature'
switch plotComponent
case 1 %an
yFreq=2*real(c);
case 2 %bn
yFreq=2*imag(c);
yFreq(1)=0;
otherwise
error('bad plotComponent combination')
end
case 'Polar'
switch plotComponent
case 1 %An
yFreq=2*abs(c);
case 2 %theta
yFreq=angle(c)*180/pi;
otherwise
error('bad plotComponent combination')
end
case 'Exponential'
switch plotComponent
case 1 %real
yFreq=real(c);
case 2 %imag
yFreq=imag(c);
case 3 %mag
yFreq=abs(c);
case 4 %angle
yFreq=angle(c)*180/pi;
otherwise
error('bad plotComponent combination')
end
end
if w<=1e-6
xFreq=zeros(size(c));
else
xFreq= linspace(0,N*w/(2*pi),N+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% Fill up output structure %%%%%%%%%%%%%%%%%%%%%%
strData.xTime = t;
strData.y0Time = y0;
strData.y1Time = y1;
strData.xFreq = xFreq;
strData.yFreq = yFreq;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ScaleIt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out=ScaleIt(sliderIn,sliderMin,sliderMean,sliderMax,outMin,outMean,outMax)
% ScaleIt provides a piecewise scaling of value SliderIn
% if SliderIn is between sliderMin and SliderMean it is scaled between outMin and outMean
% and similarly betwen Mean and Max
if sliderIn<sliderMin || sliderIn > sliderMax
error('slider in outside min/max boundaries in ScaleIt')
end
if sliderIn < sliderMean
tm = (outMean-outMin)/(sliderMean-sliderMin);
else
tm = (outMax-outMean)/(sliderMax-sliderMean);
end
tb = outMean - (sliderMean * tm);
out = sliderIn * tm + tb;