function spectrum(h, n, fo) dim = double(min(size(h)) ~= 1) + 1; switch dim case 1 spect1D(h, n, fo); case 2 spect2D(h, n, fo); end % function spect2D(h, n, fo) tagw1 = 'w1 (Hz)'; tagw2 = 'w2 (Hz)'; if isempty(fo) fo = [2 2]; tagw1 = 'w1/Pi (rad/sec)'; tagw2 = 'w2/Pi (rad/sec)'; end n1 = n(1); n2 = n(2); f1 = fo(1); f2 = fo(2); y = fftshift(fft2(h, n1, n2)); p = atan2(imag(y), real(y)); figure; freq1 = (-1:2/n1:1-2/n1)*f1/2; freq2 = (-1:2/n2:1-2/n2)*f2/2; subplot(121); contourf(freq2, freq1, abs(y)); axis ij colormap(bone) xlabel(tagw2, 'FontSize', 13); ylabel(tagw1, 'FontSize', 13); title('Magnitude Spectrum', 'FontSize', 13); subplot(122); imagesc(freq2, freq1, p/pi); axis ij colormap(bone) colorbar('EastOutside') xlabel(tagw2, 'FontSize', 13); ylabel(tagw1, 'FontSize', 13); title('Phase Spectrum', 'FontSize', 13); % function spect1D(h, n, fo) tag = 'w (Hz)'; if isempty(fo) fo = 2; tag = 'w/Pi (rad/sec)'; end w = -1:2/n:1-2/n; w = w.'; w = w*(fo/2); y = fft(h, n); y = fftshift(y); y = y(:); m = abs(y); % p = unwrap(angle(y)); p = atan2(imag(y), real(y)); figure; subplot(211); plot(w, m, 'k', 'LineWidth', 2); ylabel('Magnitude', 'FontSize', 13); axis tight title('Magnitude Spectrum', 'FontSize', 13); % phase is normalized by pi subplot(212); plot(w, p/pi, 'k', 'LineWidth', 2); ylabel('Phase/Pi', 'FontSize', 13); xlabel(tag, 'FontSize', 13); axis tight title('Phase Spectrum', 'FontSize', 13);