function reportFrequencyDistribution(this)
	import mlreportgen.report.* 
	import mlreportgen.dom.* 
	
	rpt = Report(strcat("Frequency analysis, ", inputname(1)), "docx");

	add(rpt, Heading1("Summary"));
	add(rpt, Paragraph("..."));
	
	c(:,1) = [ones(127,1); linspace(1, 0, 128)'];
	c(:,2) = [linspace(0, 1, 127)'; 1; linspace(1, 0, 127)'];
	c(:,3) = [linspace(0, 1, 128)'; ones(127,1)];

	% frequencies = [0.5 4; 4 8; 8 13; 13 30];
	frequencies = [0.5 4; 4 8];
% 	frequencies = [0.5 min(20, min([this.Data.SamplingFrequency]) / 2)];
	for iF = 1:size(frequencies, 1)
		fmin = frequencies(iF, 1);
		fn = min([this.Data.SamplingFrequency]) / 2;
		fn = frequencies(iF, 2);
		add(rpt, Heading1(sprintf("%0.1f-%0.1f Hz", fmin, fn)));
		for iFile = 1:length(this.Data)
			disp(iFile);
			add(rpt, PageBreak);
			add(rpt, Heading2(this.FileName(iFile)));
			add(rpt, Paragraph(this.MeasurementDate(iFile)));

			% dt = this.Data(iFile).Time(2) - this.Data(iFile).Time(1);
			% minscale = centfrq(wavelettype)/(fmin*dt);
			% maxscale = centfrq(wavelettype)/(fn*dt);
			% scales = logspace(log10(minscale), log10(maxscale), 100);
			% waveletCfs = [];
			% for iCell = 1:this.Data(iFile).CellNum
			% 	data = this.Data(iFile).DeltaFperF0(:,iCell,1);
			% 	waveletCfs{iCell, 1} = cwt(data, scales, wavelettype);
			% end
			% waveletFrequency = scal2frq(scales, wavelettype, dt);
			this.Data(iFile).calculateWavelet(frequencies(iF, :));

			% Wavelet summed at each frequency
			figure('units', 'centimeters', 'position', [0 0 22 16]);
			subplot(5,1,[1 3]);
			% cfs = cell2mat(cellfun(@(x) nansum(abs(x)'), waveletCfs, 'UniformOutput', false));
			% cfs = cfs ./ max(cfs(:));
			cfs = this.Data(iFile).Wavelet.CfsByFrequency;
			waveletFrequency = this.Data(iFile).Wavelet.Frequency;
			cfsg = cfs(this.Data(iFile).Glia, :);
			cfsn = cfs(this.Data(iFile).Neuron, :);
			cfsnc = cfs(~this.Data(iFile).Glia & ~this.Data(iFile).Neuron, :);
			cfs = zeros(size(cfs, 1), size(cfs, 2), 3);
			if ~isempty(cfsg)
				cfsg(isnan(cfsg)) = 0;
				cfs(1:size(cfsg, 1),:,1) = 1;
				cfs(1:size(cfsg, 1),:,2) = 1 - cfsg;
				cfs(1:size(cfsg, 1),:,3) = 1 - cfsg;
			end
			if ~isempty(cfsn)
				cfsn(isnan(cfsn)) = 0;
				cfs(size(cfsg, 1)+1:size(cfsg, 1)+size(cfsn, 1),:,1) = 1 - cfsn;
				cfs(size(cfsg, 1)+1:size(cfsg, 1)+size(cfsn, 1),:,2) = 1;
				cfs(size(cfsg, 1)+1:size(cfsg, 1)+size(cfsn, 1),:,3) = 1 - cfsn;
			end
			if ~isempty(cfsnc)
				cfsnc(isnan(cfsnc)) = 0;
				cfs(size(cfsg, 1)+size(cfsn, 1)+1:end,:,1) = 1 - cfsnc;
				cfs(size(cfsg, 1)+size(cfsn, 1)+1:end,:,2) = 1 - cfsnc;
				cfs(size(cfsg, 1)+size(cfsn, 1)+1:end,:,3) = 1;
			end
			imagesc('XData', waveletFrequency, 'YData', 1:size(cfs, 1), 'CData', cfs);
			xlim([min(waveletFrequency) max(waveletFrequency)]);
			ylim([1 this.Data(iFile).CellNum]);
			set(gca,'YDir','reverse');
			set(gca,'XScale','log');
			xticks([1 2 3 4 5 10 20 30 40 50 100]);
			title('Imaging power spectrum: astrocytes (red), neurons (green), not cell (blue)');
			% xlabel('Frequency (Hz)');
			ylabel('Cell number');

			% Wavelet summed at each frequency, summed for astrocytes or neurons
			subplot(5,1,4);
			plot(waveletFrequency, mean(cfsg), 'r');
			hold on;
			plot(waveletFrequency, mean(cfsn), 'g');
			plot(waveletFrequency, mean(cfsnc), 'b');
			set(gca, 'XScale', 'log');
			title('Total wavelet coefficient for all cells');
			xlabel('Frequency (Hz)');

			% Ephys wavelet summed at each frequency
			if ~isempty(this.Data(iFile).Ephys)
				wavelettype = 'cmor1-2';
				dtE = (this.Data(iFile).EphysTime(2) - this.Data(iFile).EphysTime(1));
				minscale = centfrq(wavelettype)/(max(fmin,0.5)*dtE);
				maxscale = centfrq(wavelettype)/(fn*dtE);
				scales = logspace(log10(minscale), log10(maxscale), 100);
				data = this.Data(iFile).Ephys;
				ephysWaveletCfs = cwt(data, scales, wavelettype);
				ephysWaveletFrequency = scal2frq(scales, wavelettype, dtE);

				subplot(3,1,3);
				cfs = nansum(abs(ephysWaveletCfs)');
				plot(ephysWaveletFrequency, cfs);
				xlim([fmin fn]);
				set(gca,'XScale','log');
				% xticks([0.5 1 2 3 4 5 10 20 30 40 50 100]);
				title('Ephys power spectrum');
				xlabel('Frequency (Hz)');
				ylabel('Sum power');
			end

			spaceplots;
			addFigureToReport(rpt, 16, 20);
			close(gcf);

			% Phase lock
			figure('units', 'centimeters', 'position', [0 0 10 16]);
			this.Data(iFile).calculatePhaselock(frequencies(iF, :));
			
			plv = this.Data(iFile).PhaseLock;
			subplot(5,2,1);
			imagesc(plv);
			clim([0 1]);
			title("PLV, all cells");

			plvg = this.Data(iFile).PhaseLock(this.Data(iFile).Glia, this.Data(iFile).Glia);
			subplot(5,2,2);
			imagesc(plvg);
			clim([0 1]);
			title("PLV, astrocytes");

			plvn = this.Data(iFile).PhaseLock(this.Data(iFile).Neuron, this.Data(iFile).Neuron);
			subplot(5,2,3);
			imagesc(plvn);
			clim([0 1]);
			title("PLV, neurons");
			
			plvnc = this.Data(iFile).PhaseLock(~this.Data(iFile).Glia & ~this.Data(iFile).Neuron, ~this.Data(iFile).Glia & ~this.Data(iFile).Neuron);
			subplot(5,2,4);
			imagesc(plvnc);
			clim([0 1]);
			title("PLV, not cells");

			plvgn = this.Data(iFile).PhaseLock(this.Data(iFile).Glia, this.Data(iFile).Neuron);
			subplot(5,2,5);
			imagesc(plvgn);
			clim([0 1]);
			title("PLV, astrocyte to neurons");

			subplot(5,2,[7 8]);
			plvg = tril(plvg);
			plvg = plvg(plvg>0);
			plvn = tril(plvn);
			plvn = plvn(plvn>0);
			plvnc = tril(plvnc);
			plvnc = plvnc(plvnc>0);
			plvgn = tril(plvgn);
			plvgn = plvgn(plvgn>0);
			plvall = [plvg; plvn; plvnc; plvgn];
			celltypeall = [repmat("astrocytes", length(plvg), 1); repmat("neurons", length(plvn), 1); repmat("not cell", length(plvnc), 1); repmat("astrocyte to neuron", length(plvgn), 1)];
			swarmchart(categorical(celltypeall), plvall, 8, 'filled', 'XJitter', 'rand', 'MarkerFaceAlpha', 0.5, 'MarkerEdgeAlpha', 0.5);
			ylim([0 1]);
			
			subplot(5,2,[9 10]);
			histogram(plvg, 'DisplayStyle', 'stair', 'EdgeColor', 'red');
			hold on;
			histogram(plvn, 'DisplayStyle', 'stair', 'EdgeColor', 'green');
			histogram(plvnc, 'DisplayStyle', 'stair', 'EdgeColor', 'blue');
			histogram(plvgn, 'DisplayStyle', 'stair', 'EdgeColor', 'black');

			% spaceplots;
			addFigureToReport(rpt, 16, 22);
			close(gcf);
		end
	end
		

	close(rpt);
end