Subversion Repositories wbpwmaudio

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Filename: 	showspectrogram.m
%%
%% Project:	A Wishbone Controlled PWM (audio) controller
%%
%% Purpose:	To generate a spectrogram image which can then be used to
%%		evaluate the wavfp.dbl file produced by the pdmdemo executable
%%	in this directory.
%%
%%	This script has been tested with Octave, and it is known to work with
%%	Octave.  While it may work within Matlab as well, no representation is
%%	being made to that effect.
%%
%% Creator:	Dan Gisselquist, Ph.D.
%%		Gisselquist Technology, LLC
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Copyright (C) 2017, Gisselquist Technology, LLC
%%
%% This program is free software (firmware): you can redistribute it and/or
%% modify it under the terms of the GNU General Public License as published
%% by the Free Software Foundation, either version 3 of the License, or (at
%% your option) any later version.
%%
%% This program is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTIBILITY or
%% FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
%% for more details.
%%
%% You should have received a copy of the GNU General Public License along
%% with this program.  (It's in the $(ROOT)/doc directory.  Run make with no
%% target there if the PDF file isn't present.)  If not, see
%% <http://www.gnu.org/licenses/> for a copy.
%%
%% License:	GPL, v3, as defined and found on www.gnu.org,
%%		http://www.gnu.org/licenses/gpl.html
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%
function [img] = showspectrogram(fname)
	% Read the waveform into memory
	fid = fopen(fname,'r');
	dat = fread(fid, inf, 'double');
	fclose(fid);
 
	% You can plot it at this point if you would like
	% plot(dat);
 
	sample_rate   = 44.1e3;
	window_length = 512; % Size of the FFT
	Nc = 256; % Number of colors
 
	h = hanning(window_length); % About 1 ms
	rows = (length(dat)-length(h)) / length(h) * 2
	img = zeros(length(h)/2, rows);
	size(img)
 
	for n=1:rows
		cut = dat((1:length(h)) + (n-1)*length(h)/2);
		cutf = abs(fft(cut.*h)).^2;
		img(:,n) = cutf(1:(length(h)/2));
	end
 
	mintime = 0;
	maxtime = rows / sample_rate * window_length / 2;
	minfreq = 0;
	maxfreq = 0.5 * sample_rate;
 
	% Normalize the image so that the maximum value is one.
	img = img ./ max(max(img));
 
	%
	% Turn out spectrogram image into dB
	%
	% Adding 1e-8 is useful for artificially forcing the floor of the
	% image to be at a particular value, as well as for keeping the log
	% from reporting values too negative to plot successfully.
	img = 10 * log(img + 1e-8)/log(10);
 
	img = (img + 80)/80;
	colormap(mymap(256));
	image([mintime, maxtime], [minfreq, maxfreq/1e3], img * Nc);
	ylabel('Frequency (kHz)');
	xlabel('Time (s)');
 
	% Trim this image output to the relevant portion of the output
	axis([0, 6.2, 0, maxfreq/1e3]);
 
% For the plot, comment out the img ./ max(max(img)) line, as well as the
% img = (img+80)/80 line.  Then you can call this as:
% img    = showspectrogram('pdm8k-weak.dbl');
% imgpwm = showspectrogram('pwm8k-weak.dbl');
%
% freq = (0:(window_length-1))/(window_length)*(sample_rate/2/1e3)
% plot(freq,img(:,1400),'b;PDM;',freq,imgpwm(:,1400),'r;PWM;');
% grid on; xlabel('Frequency (kHz)'); ylabel('Volume (dB)');
% axis([0,22,-70,50]);