spectrumexpofn.m

非常好的数字处理教程

%
function spectrumexpofn(action)

    if nargin < 1
        action='init';
    end

	% Fix limits

	name = mfilename;
    figname = [name(1:end-2) '_fig'];
    f=findobj('Tag',figname);
    handles = get(f,'UserData');

    switch action
        case 'help'
			display_help(figname);
		case 'init'
            setdefaults;
			reset(handles.sin1);
			reset(handles.sinsum);
			reset(handles.freqdomain);
			reset(handles.phasdomain);
			update(handles);
		case 'update'
			update(handles);
		case 'play'
			[signal, Fs] = calc_signal(handles);
			signal = normalize(signal);
			audiodata.data = signal';
			audiodata.Fs = Fs;
			play_audiodata(audiodata, handles.play);
        case {'fourier','sonogram','alias'}
			[signal, Fs] = calc_signal(handles);
			signal = normalize(signal);
			audiodata.data = signal';
			audiodata.Fs = Fs;
			if (max(signal) > 1.0)
				signal = normalize(signal);
			end
			switch action
				case 'fourier'
					fourierexpogui(audiodata);
				case 'sonogram'
					sonoexpogui(audiodata);
                case 'alias'
                    aliasexpogui(audiodata);
			end
			
	end
	set(f,'UserData',handles);


% --------------------------------------------------------------------
function update(handles)
	axes(handles.sin1);
	cla;
	axes(handles.sinsum);
	cla;
	axes(handles.freqdomain);
	cla;
	axes(handles.phasdomain);
	cla;

	Fs = 44100;
	duration = str2double(get(handles.durationtext,'String'));
	t = [0:1/Fs:duration];
	freq1 = str2double(get(handles.freq1,'String'));
	freq2 = str2double(get(handles.freq2,'String'));
	amp1 = str2double(get(handles.amp1,'String'));
	amp2 = str2double(get(handles.amp2,'String'));
	phase1 = str2double(get(handles.phase1,'String'))*pi/180;
	phase2 = str2double(get(handles.phase2,'String'))*pi/180;
	offset = str2double(get(handles.offset,'String'));
	sine1 = offset + amp1*cos(2*pi*freq1*t + phase1);
	sine2 = amp2*cos(2*pi*freq2*t + phase2);

	axes(handles.sin1);
	if (get(handles.plot1box,'Value') & get(handles.plot2box,'Value'))
		plot(t,sine1,'b-',t,sine2,'k--');
		legend('sine1','sine2');
		hold on;
	elseif (get(handles.plot1box,'Value'))
		plot(t,sine1);
		hold on;
		legend('sine1');
	elseif (get(handles.plot2box,'Value')) 
		plot(t,sine2);
		hold on;
		legend('sine2');
	end

	popup_sel_index = get(handles.popupmenu1, 'Value');
	switch popup_sel_index
		case 1
			% Sum
			signal = sine1 + sine2;
			legendtext = 'sine1 + sine2';
		case 2
			% Product
			signal = sine1.*sine2;
			legendtext = 'sine1*sine2';
	end

	axes(handles.sinsum);
	plot(t,signal);
	legend(legendtext);
	xlabel('Time (s)');
	ylabel('Amplitude');
    
	switch popup_sel_index
		case 1
			% Sum
			axes(handles.freqdomain);
			set(handles.freqdomain,'XLimMode','auto');
			set(handles.freqdomain,'YLimMode','auto');
			if get(handles.plot1box,'Value')
				stem([-freq1 freq1],[amp1/2 amp1/2]);
				if offset > 0
					hold on;
					stem(0,offset);
				end
			end
			hold on;
			if get(handles.plot2box,'Value')
				stem([-freq2 freq2],[amp2/2 amp2/2]);
			end
			maxfreq = max([freq1 freq2]);
			maxamp = sum([amp1 amp2]);
			if freq1 == freq2
				cla;
				amp = (amp1*cos(phase1) + amp2*cos(phase2))^2 + ...
						(amp1*sin(phase1) + amp2*sin(phase2))^2;
				amp = sqrt(amp);
				stem([-freq1 freq1],[amp/2 amp/2]);
				if offset > 0
					hold on;
					stem(0,offset);
				end
				maxamp = max([amp amp2]);
			end
			xticks = get(handles.freqdomain,'XTick');
			yticks = get(handles.freqdomain,'YTick');
			xlim([1.1*xticks(1) 1.1*xticks(end)]);
			ylim([0 1.1*yticks(end)]);

			axes(handles.phasdomain);
			set(handles.phasdomain,'XLimMode','auto');
			set(handles.phasdomain,'YLimMode','auto');
			if get(handles.plot1box,'Value')
				stem([-freq1 freq1],[-phase1*180/pi phase1*180/pi]);
			end
			hold on;
			if get(handles.plot2box,'Value')
				stem([-freq2 freq2],[-phase2*180/pi phase2*180/pi]);
			end
			if freq1 == freq2
				cla;
				phase = (amp1*sin(phase1) + amp2*sin(phase2))/...
						(amp1*cos(phase1) + amp2*cos(phase2));
				phase = atan(phase);
				stem([-freq1 freq1],[-phase*180/pi phase*180/pi]);
			end
			xlim([1.1*xticks(1) 1.1*xticks(end)]);
		case 2
			% Product
			set(handles.freqdomain,'XLimMode','auto');
			set(handles.freqdomain,'YLimMode','auto');
			axes(handles.freqdomain);
			amp = amp1*amp2;
			freqs = [freq1+freq2 freq1-freq2];
			freqs = [-freqs freqs];
			stem(freqs, [amp/4 amp/4 amp/4 amp/4]);
			hold on;
			if offset > 0
				stem([-freq2 freq2], [offset*amp2/2 offset*amp2/2]);
			end
			xticks = get(handles.freqdomain,'XTick');
			yticks = get(handles.freqdomain,'YTick');
			xlim([1.1*xticks(1) 1.1*xticks(end)]);
			ylim([0 1.1*yticks(end)]);

			axes(handles.phasdomain);
			set(handles.phasdomain,'XLimMode','auto');
			set(handles.phasdomain,'YLimMode','auto');
			freqs = [freq1+freq2 freq1-freq2];
			freqs = [-freqs freqs];
			phas = [phase1+phase2 phase1-phase2];
			phas = [-phas.*180/pi phas.*180/pi];
			stem(freqs, phas);
			hold on;
			if offset > 0
				stem([-freq2 freq2], [-phase2.*180/pi phase2.*180/pi]);
			end 
			xlim([1.1*xticks(1) 1.1*xticks(end)]);
	end
	axes(handles.freqdomain);
	ylabel('Magnitude');

	axes(handles.phasdomain);
	xlabel('Frequency (Hz)');
	ylabel('Phase (deg)');

	if (get(handles.gridbox,'Value'))
		axes(handles.sin1);
		grid on;
		axes(handles.sinsum);
		grid on;
		axes(handles.phasdomain);
		grid on;
		axes(handles.freqdomain);
		grid on;
	else
		axes(handles.sin1);
		grid off;
		axes(handles.sinsum);
		grid off;
		axes(handles.phasdomain);
		grid off;
		axes(handles.freqdomain);
		grid off;
	end

	% Set plot properties
	set(handles.sin1,'XTickLabel',['']);
	set(handles.freqdomain,'XTickLabel',['']);
	set(handles.sin1,'YLim',[-1.5 1.5]);
	set(handles.sinsum,'YLim',[-2 2]);
	set(handles.sin1,'XLim',[0 duration]);
	set(handles.sinsum,'XLim',[0 duration]);

% --------------------------------------------------------------------
function [signal,Fs] = calc_signal(handles)
	Fs = 44100;
	freq1 = str2double(get(handles.freq1,'String'));
	freq2 = str2double(get(handles.freq2,'String'));
	amp1 = str2double(get(handles.amp1,'String'));
	amp2 = str2double(get(handles.amp2,'String'));
	phase1 = str2double(get(handles.phase1,'String'))*pi/180;
	phase2 = str2double(get(handles.phase2,'String'))*pi/180;
	offset = str2double(get(handles.offset,'String'));
	t = [0:1/Fs:2];
	sine1 = offset + amp1*sin(2*pi*freq1*t + phase1);
	sine2 = amp2*sin(2*pi*freq2*t + phase2);
	popup_sel_index = get(handles.popupmenu1, 'Value');
	switch popup_sel_index
		case 1
			% Sum
			signal = sine1 + sine2;
		case 2
			% Product
			signal = sine1.*sine2;
	end