www.gusucode.com > sigdemos 工具箱matlab源码程序 > sigdemos/DFTExample.m
function DFTExample(action)
%DFTExample Interactive DFT of a signal
% This example shows how to use some basic signal processing concepts,
% such as sampling, aliasing, and windowing.
%
% You are seeing discrete samples of a periodic waveform (the upper plot) and
% the absolute value of its discrete Fourier transform (DFT), obtained using
% a fast Fourier transform (FFT) algorithm (the lower plot).
%
% In the lower plot, frequencies from 0 to 100 Hertz are displayed. The DFT
% at negative frequencies is a mirror image of the DFT at positive frequen-
% cies. The sampling rate is 200 Hertz, which means the "Nyquist frequency"
% is 100 Hertz. The DFT at frequencies above the Nyquist frequency is the
% same as the DFT at lower (negative) frequencies.
%
% Click and drag a point on the waveform shown in the upper plot to move that
% point to a new location. This sets a new fundamental frequency and
% amplitude.
%
% The "Signal" pop-up menu lets you change the shape of the waveform.
%
% The "Window" menu lets you select a window function. This window is multi-
% plied by the time waveform prior to taking the DFT.
%
% The fundamental frequency of the waveform is given in the text box. You can
% change the fundamental frequency by clicking in the text box, editing the
% number there, and pressing RETURN. You can also change the fundamental by
% clicking and dragging on the waveform.
% Author: T. Krauss
% 11/3/92, updated 2/9/93
% Adapted for Expo: dlc, 7-93
% Copyright 1988-2012 The MathWorks, Inc.
% possible actions:
% 'start'
% 'down'
% 'move'
% 'up'
% 'redraw'
% 'done'
% 'setfreq'
% 'setwindow'
% 'showwind'
if nargin<1,
action='start';
end;
global SIGDEMO1_DAT
global ADDIT_DAT % this is for WINDOW pop-up button
if strcmp(action,'start'),
%====================================
% Graphics initialization
oldFigNumber = watchon;
figNumber = figure;
set(figNumber, ...
'NumberTitle','off', ...
'Name','Discrete Fourier Transform', ...
'Units','normalized');
%====================================
% Information for all buttons
labelColor=192/255*[1 1 1];
top=0.95;
bottom=0.05;
yInitLabelPos=0.90;
left = 0.78;
labelWid=0.18;
labelHt=0.05;
btnWid = 0.18;
btnHt=0.07;
% Spacing between the label and the button for the same command
btnOffset=0.003;
% Spacing between the button and the next command's label
spacing=0.05;
%====================================
% The CONSOLE frame
frmBorder=0.02;
yPos=0.05-frmBorder;
frmPos=[left-frmBorder yPos btnWid+2*frmBorder 0.9+2*frmBorder];
h=uicontrol( ...
'Style','frame', ...
'Units','normalized', ...
'Position',frmPos);
%====================================
% The SIGNAL command popup button
btnNumber=1;
yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing);
labelStr=' Signal';
% Generic label information
labelPos=[left yLabelPos-labelHt labelWid labelHt];
uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'HorizontalAlignment','left', ...
'String',labelStr);
btnPos=[left yLabelPos-labelHt-btnHt-btnOffset btnWid btnHt];
popup=uicontrol('Style','Popup','String','sine|square|sawtooth',...
'Units','normalized',...
'Position', btnPos, ...
'Callback','DFTExample(''redraw'')');
%====================================
% The WINDOW command popup button
btnNumber=2;
yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing);
labelStr=' Window';
popupStr= ' rectangle| triangular| hanning| hamming| chebyshev| kaiser';
callbackStr= 'DFTExample(''setwindow'')';
% Generic label information
labelPos=[left yLabelPos-labelHt labelWid labelHt];
uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position',labelPos, ...
'HorizontalAlignment','left', ...
'String',labelStr);
% Generic popup button information
btnPos=[left yLabelPos-labelHt-btnHt-btnOffset btnWid btnHt];
winHndl = uicontrol( ...
'Style','popup', ...
'Units','normalized', ...
'Position',btnPos, ...
'String',popupStr, ...
'Callback',callbackStr);
%====================================
% The FUNDAMENTAL editable text box
btnNumber=3;
yLabelPos=top-(btnNumber-1)*(btnHt+labelHt+spacing);
labelPos=[left yLabelPos-labelHt labelWid labelHt];
freq_text = uicontrol( ...
'Style','text', ...
'Units','normalized', ...
'Position', labelPos, ...
'String','Fundamental');
btnPos=[left+0.02 yLabelPos-labelHt-btnHt-btnOffset ...
0.5*btnWid+frmBorder btnHt];
freq_field = uicontrol( ...
'Style','edit', ...
'Units','normalized', ...
'Position', btnPos, ...
'BackgroundColor','w',...
'String',' 5',...
'Callback','DFTExample(''setfreq''); DFTExample(''redraw'')');
%====================================
% The INFO button
uicontrol( ...
'Style','push', ...
'Units','normalized', ...
'Position',[left bottom+(2*labelHt)+spacing btnWid 2*labelHt], ...
'String','Info', ...
'Callback','DFTExample(''info'')');
%========================================
% The CLOSE button
done_button=uicontrol('Style','Pushbutton', ...
'Units','normalized',...
'Position',[left bottom btnWid 2*labelHt], ...
'Callback','DFTExample(''done'')','String','Close');
%====================================
% Create initial signal
N=201; % number of samples
% N=21;
amp=0.5;
freq=5; % hertz
t0=0; % seconds
t1=1; % seconds
t=linspace(t0,t1,N)';
T=(t1-t0)/N; % sampling rate in seconds
M=256; % length of fft
window=ones(N,1); % use window to lower side-lobes in the freq. domain
% (makes peaks wider)
min_dB = -40; % freq. domain lower axis limit
% create axes for time domain and frequency domain plot
ax_freq=axes('Position',[.12 .14 .6 .3],'XLim',...
[0 1/(2*T)],'YLim',[min_dB 50]);
val = get(popup,'Value');
% time domain
f=local_dataSet(val,amp,t,freq);
% frequency domain
[FF,w]=local_fftpsd(f,window,M,min_dB);
freq_line=plot(w/2/pi/T,FF);
axis([0 1/(2*T) min_dB 50]);
grid on;
ylabel('Magnitude (dB)');
xlabel('Frequency (Hertz)');
ax_time=axes('Position',[.12 .58 .6 .3],'XLim',[t0 t1],'YLim',[-1 1]);
time_line=plot(t,f);
axis([t0 t1 -1 1]);
% (set to xor mode to prevent re-rendering, that is, for speed)
grid on;
ylabel('Waveform');
xlabel('Time (Seconds)');
text(0.12, 1.55,' Click and drag waveform to change');
text(0.12, 1.3,'fundamental frequency and amplitude');
set(time_line,'ButtonDownFcn','DFTExample(''down'')');
SIGDEMO1_DAT = {freq; amp; N; M; min_dB; 0; 0; ...
time_line; freq_line; freq_field; popup; -1; gcf; t(:); window(:)};
ADDIT_DAT = winHndl;
watchoff(oldFigNumber);
elseif strcmp(action,'down'),
% assumes that a line was clicked
time_line=SIGDEMO1_DAT{8};
axes(get(time_line,'Parent')); % set the right axes
% Obtain coordinates of mouse click location in axes units
pt=get(gca,'CurrentPoint');
x=pt(1,1);
y=pt(1,2);
% find closest vertex of line to mouse click location (call it fixed_x, fixed_y)
line_x=get(time_line,'XData');
line_y=get(time_line,'YData');
units_str = get(gca,'Units'); % save normalized state
set(gca,'Units','pixels'); % distance must be in pixels
p=get(gca,'Position');
xa=get(gca,'XLim');
ya=get(gca,'YLim');
dist=((line_x-x)*p(3)/(xa(2)-xa(1))).^2 + ...
((line_y-y)*p(4)/(ya(2)-ya(1))).^2;
[temp,i]=min(dist);
fixed_x=line_x(i);
fixed_y=line_y(i);
set(time_line,'LineStyle','--');
SIGDEMO1_DAT{6}=fixed_x;
SIGDEMO1_DAT{7}=fixed_y;
set(gca,'Units',units_str );
set(gcf,'WindowButtonMotionFcn', 'DFTExample(''move'')');
set(gcf,'WindowButtonUpFcn', 'DFTExample(''up'')');
% set(gcf,'userdata',u);
elseif strcmp(action,'move'),
% u = get(gcf,'userdata');
freq=SIGDEMO1_DAT{1};
amp=SIGDEMO1_DAT{2};
N=SIGDEMO1_DAT{3};
M=SIGDEMO1_DAT{4};
min_dB=SIGDEMO1_DAT{5};
fixed_x=SIGDEMO1_DAT{6};
fixed_y=SIGDEMO1_DAT{7};
time_line=SIGDEMO1_DAT{8};
freq_line=SIGDEMO1_DAT{9};
freq_field=SIGDEMO1_DAT{10};
popup=SIGDEMO1_DAT{11};
t=SIGDEMO1_DAT{14};
window=SIGDEMO1_DAT{15};
% Avoid divide by zero warnings
if fixed_y==0,
fixed_y=eps;
end
pt=get(gca,'CurrentPoint');
x=pt(1,1);
y=pt(1,2);
% Avoid divide by zero warnings
if x==0,
x=eps;
end
amp1=y/fixed_y*amp;
if (abs(amp1)>1.0),
amp1=1.0*sign(amp1);
end;
if (abs(amp1)<0.05),
amp1=0.05*sign(amp1);
end;
if (amp1 == 0),
amp1=0.05;
end;
freq1=fixed_x/x*freq;
val = get(popup,'Value');
% time domain
f=local_dataSet(val,amp1,t,freq1);
set(time_line,'YData',f);
% frequency domain
[FF,w]=local_fftpsd(f,window,M,min_dB);
set(freq_line,'YData',FF);
set(freq_field,'String',num2str(freq1));
elseif strcmp(action,'up'),
pt=get(gca,'CurrentPoint');
x=pt(1,1);
y=pt(1,2);
set(gcf,'WindowButtonMotionFcn','');
set(gcf,'WindowButtonUpFcn','');
% u=get(gcf,'userdata');
freq=SIGDEMO1_DAT{1};
amp=SIGDEMO1_DAT{2};
fixed_x=SIGDEMO1_DAT{6};
fixed_y=SIGDEMO1_DAT{7};
if fixed_y==0,
fixed_y=eps;
end
amp1=y/fixed_y*amp;
if (abs(amp1)>1.0),
amp1=1.0*sign(amp1);
end;
if (abs(amp1)<0.05),
amp1=0.05*sign(amp1);
end;
freq1=fixed_x/x*freq;
set(SIGDEMO1_DAT{8},'LineStyle','-');
SIGDEMO1_DAT{1}=freq1; % set amplitude and frequency
SIGDEMO1_DAT{2}=amp1;
% set(gcf,'userdata',u);
DFTExample('redraw');
elseif strcmp(action,'done'),
% close the figure window that is showing the window function:
% u = get(gcf,'userdata');
if (SIGDEMO1_DAT{12}~=-1),
close(SIGDEMO1_DAT{12});
end;
close(gcf);
clear global SIGDEMO1_DAT
clear global ADDIT_DAT
elseif strcmp(action,'redraw'),
% recomputes time and frequency waveforms and updates display
% u = get(gcf,'userdata');
freq=SIGDEMO1_DAT{1};
amp=SIGDEMO1_DAT{2};
N=SIGDEMO1_DAT{3};
M=SIGDEMO1_DAT{4};
min_dB=SIGDEMO1_DAT{5};
time_line=SIGDEMO1_DAT{8};
freq_line=SIGDEMO1_DAT{9};
freq_field=SIGDEMO1_DAT{10};
popup=SIGDEMO1_DAT{11};
t=SIGDEMO1_DAT{14};
window=SIGDEMO1_DAT{15};
val = get(popup,'Value');
% time domain
f=local_dataSet(val,amp,t,freq);
set(time_line,'YData',f);
% frequency domain
[FF,w]=local_fftpsd(f,window,M,min_dB);
set(freq_line,'YData',FF);
set(freq_field,'String',num2str(freq));
drawnow;
elseif strcmp(action,'setwindow'),
% u = get(gcf,'userdata');
winHndl = ADDIT_DAT;
in1 = get(winHndl,'Value');
in2 = 30;
N=SIGDEMO1_DAT{3};
if (in1==1),
window = boxcar(N);
elseif (in1==2),
window = triang(N);
elseif (in1==3),
window = hanning(N);
elseif (in1==4),
window = hamming(N);
elseif (in1==5),
window = chebwin(N,30);
elseif (in1==6),
window = kaiser(N,4);
end;
SIGDEMO1_DAT{15}=window;
% set(gcf,'userdata',u);
DFTExample('redraw');
if (SIGDEMO1_DAT{12}~=-1),
DFTExample('showwind');
end;
elseif strcmp(action,'showwind'),
% u=get(gcf,'userdata');
oldfig=gcf;
N=SIGDEMO1_DAT{3};
t=SIGDEMO1_DAT{14};
window=SIGDEMO1_DAT{15};
if (SIGDEMO1_DAT{12}==-1),
SIGDEMO1_DAT{12}=figure;
axes('Position',[.15 .62 .8 .3]);
line1=plot(t,window);
title('Window function');
xlabel('time (seconds)');
grid on;
ylabel('Window');
axes('Position',[.15 .2 .8 .3]);
W=fft(window,1024);
line2=plot((0:(1/1024):(.5-(1/1024)))*N,20*log10(abs(W(1:512))));
set(gca,'XLim',[0 N/2]);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;
windclose=uicontrol('Style','Pushbutton','Units','Normalized',...
'Position',[.85 .02 .12 .08],...
'Callback',['THEHAND=get(gcf,''UserData'');'...
'close; figure(THEHAND(1)); global SIGDEMO1_DAT, '...
'SIGDEMO1_DAT{12}=-1; clear THEHAND SIGDEMO1_DAT'],...
'String','Close');
set(gcf,'UserData',[oldfig line1 line2]);
% set(oldfig,'userdata',u);
else
figure(SIGDEMO1_DAT{12});
lines=get(gcf,'UserData');
set(lines(2),'YData',window);
W=fft(window,1024);
set(lines(3),'YData',20*log10(abs(W(1:512))));
drawnow
end
elseif strcmp(action,'setfreq'),
x = str2double(get(SIGDEMO1_DAT{10},'String'));
if isnan(x), % handle the non-numeric case
set(SIGDEMO1_DAT{10},'String',num2str(SIGDEMO1_DAT{1}));
else
SIGDEMO1_DAT{1}=x;
end;
elseif strcmp(action,'info'),
ttlStr = 'Discrete Fourier Transform';
hlpStr= ...
['You are seeing discrete samples of a periodic waveform (the upper plot) and'
'the absolute value of its discrete Fourier transform (DFT), obtained using '
'a fast Fourier transform (FFT) algorithm (the lower plot). '
' '
'In the lower plot, frequencies from 0 to 100 Hertz are displayed. The DFT '
'at negative frequencies is a mirror image of the DFT at positive frequen- '
'cies. The sampling rate is 200 Hertz, which means the "Nyquist frequency" '
'is 100 Hertz. The DFT at frequencies above the Nyquist frequency is the '
'same as the DFT at lower (negative) frequencies. '
' '
'Click and drag a point on the waveform shown in the upper plot to move that'
'point to a new location. This sets a new fundamental frequency and '
'amplitude. '
' '
'The "Signal" pop-up menu lets you change the shape of the waveform. '
' '
'The "Window" menu lets you select a window function. This window is multi- '
'plied by the time waveform prior to taking the DFT. '
' '
'The fundamental frequency of the waveform is given in the text box. You can'
'change the fundamental frequency by clicking in the text box, editing the '
'number there, and pressing RETURN. You can also change the fundamental by '
'clicking and dragging on the waveform. '];
helpwin(hlpStr, ttlStr);
end
%-------------------------------------------------------------
function [FF,w] = local_fftpsd(f,window,M,min_dB)
% LOCAL_FFTPSD Calculates the Power Spectral Density estimate
% via the FFT.
%
F=fft(window.*f,2*M);
F=F(1:M);
w=(0:M-1)*pi/M;
indxs = find(F==0);
F(indxs) = eps; % Avoid warning of divide by zero
FF=20*log10(abs(F));
ind=find(FF<min_dB);
FF(ind)=NaN*ones(size(ind)); % put NaN's in where
% min_dB shows up - this is to work around no clipping in xor mode
%-------------------------------------------------------------
function f = local_dataSet(val,amp,t,freq)
% LOCAL_DATASET - Calculates the test data set based on the signal
% the user selected in the popupmenu.
if (val == 1),
f=amp*sin(freq*t*2*pi);
elseif (val == 2), % square wave
tt=freq*t*2*pi;
tmp=rem(tt,2*pi);
f=amp*(2*rem((tt<0)+(tmp>pi | tmp<-pi)+1,2)-1);
elseif (val == 3), % sawtooth
tt=freq*t*2*pi;
f=amp*((tt < 0) + rem(tt,2*pi)/2/pi - .5)*2;
end;