function spectrum2(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 = 'f1 (Hz)'; tagw2 = 'f2 (Hz)'; if isempty(fo) fo = [1 1]; 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 = (-0.5:1/n1:0.5-1/n1)*f1; freq2 = (-0.5:1/n2:0.5-1/n2)*f2; 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 = 'f (Hz)'; if isempty(fo) fo = 1; end w = -0.5:1/n:0.5-1/n; w = w.'; w = w*fo; 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);