MATLAB | πDay快乐，pi居然有这么多炫酷好玩的可视化

Posted 2023-03-15 slandarer

πDay快乐！

3.14 π Day 到了，本期给出一些非常炫酷的pi的可化并给出完整 MATLAB 代码。

首先,我们要获得pi的前n位小数，需要编写如下代码(为了防止不精准会多取几位尾数，用的时候再进行一次取前多少位的操作)：

function Pi=getPi(n)
if nargin<1,n=3;end
Pi=char(vpa(sym(pi),n+10));
Pi=abs(Pi)-48;
Pi=Pi(3:end-4);
end 

有了这个获取pi小数位的函数，我们的可视化之旅正式启动_{从简单到复杂一步步来}(请尽量使用比较新的版本运行代码，至少17b之后版本)

---

1 饼图

就是统计一下前1500位小数各个数字所占比例：

% 获取pi前1500位小数
Pi=getPi(1500);
% 统计各个数字出现次数
numNum=find([diff(sort(Pi)),1]);
numNum=[numNum(1),diff(numNum)];
% 配色列表
CM=[20,164,199;43,187,170;53,165,81;189,190,28;248,167,22;
    232,74,27;244,57,99;240,118,177;168,109,195;78,125,187]./255;
% 绘图并修饰
pieHdl=pie(numNum);
set(gcf,'Color',[1,1,1],'Position',[200,100,620,620]);
for i=1:2:20
    pieHdl(i).EdgeColor=[1,1,1];
    pieHdl(i).LineWidth=1;
    pieHdl(i).FaceColor=CM((i+1)/2,:);
end
for i=2:2:20
    pieHdl(i).Color=CM(i/2,:);
    pieHdl(i).FontWeight='bold';
    pieHdl(i).FontSize=14;
end
% 绘制图例并修饰
lgdHdl=legend(num2cell('0123456789'));
lgdHdl.FontWeight='bold';
lgdHdl.FontSize=11;
lgdHdl.TextColor=[.5,.5,.5];
lgdHdl.Location='southoutside';
lgdHdl.Box='off';
lgdHdl.NumColumns=10;
lgdHdl.ItemTokenSize=[20,15];
title("VISUALIZING  \\pi 'Pi' Chart | 1500 digits",'FontSize',18,...
    'FontName','Times New Roman','Color',[.5,.5,.5])

---

2 折线图

计算每种数字所占比例的变化：

% 获取pi前1500位小数
Pi=getPi(1500); 
% 计算比例变化
Ratio=cumsum(Pi==(0:9)',2);
Ratio=Ratio./sum(Ratio);
D=1:length(Ratio);
% 配色列表
CM=[20,164,199;43,187,170;53,165,81;189,190,28;248,167,22;
    232,74,27;244,57,99;240,118,177;168,109,195;78,125,187]./255;
hold on
% 循环绘图
for i=1:10
    plot(D(20:end),Ratio(i,20:end),'Color',[CM(i,:),.6],'LineWidth',1.8)
end
% 坐标区域修饰
ax=gca;box on;grid on
ax.YLim=[0,.2];
ax.YTick=0:.05:.2;
ax.XTick=0:200:1400;
ax.YTickLabel='0%','5%','10%','15%','20%';
ax.XMinorTick='on';
ax.YMinorTick='on';
ax.LineWidth=.8;
ax.GridLineStyle='-.';
ax.FontName='Cambria';
ax.FontSize=11;
ax.XLabel.String='Decimals';
ax.YLabel.String='Proportion';
ax.XLabel.FontSize=13;
ax.YLabel.FontSize=13;
% 绘制图例并修饰
lgdHdl=legend(num2cell('0123456789'));
lgdHdl.NumColumns=5;
lgdHdl.FontWeight='bold';
lgdHdl.FontSize=11;
lgdHdl.TextColor=[.5,.5,.5];

---

3 堆叠填充图

计算每种数字所占比例的变化：

% 获取pi前500位小数
Pi=getPi(500); 
% 计算比例变化
Ratio=cumsum(Pi==(0:9)',2);
Ratio=Ratio./sum(Ratio);
% 配色列表
CM=[231,98,84;239,138,71;247,170,88;255,208,111;255,230,183;
    170,220,224;114,188,213;82,143,173;55,103,149;30,70,110]./255;
% 绘制堆叠面积图
hold on
areaHdl=area(Ratio');
for i=1:10
    areaHdl(i).FaceColor=CM(i,:);
    areaHdl(i).FaceAlpha=.9;
end
% 图窗和坐标区域修饰
set(gcf,'Position',[200,100,720,420]);
ax=gca;
ax.YLim=[0,1];
ax.XMinorTick='on';
ax.YMinorTick='on';
ax.LineWidth=.8;
ax.FontName='Cambria';
ax.FontSize=11;
ax.TickDir='out';
ax.XLabel.String='Decimals';
ax.YLabel.String='Proportion';
ax.XLabel.FontSize=13;
ax.YLabel.FontSize=13;
ax.Title.String='Area Chart of Proportion — 500 digits';
ax.Title.FontSize=14;
% 绘制图例并修饰
lgdHdl=legend(num2cell('0123456789'));
lgdHdl.NumColumns=5;
lgdHdl.FontSize=11;
lgdHdl.Location='southeast';

---

4 连接堆叠柱状图

还是计算每种数字所占比例的变化：

% 获取pi前100位小数
Pi=getPi(100); 
% 计算比例变化
Ratio=cumsum(Pi==(0:9)',2);
Ratio=Ratio./sum(Ratio);

X=Ratio(:,10:10:80)';
barHdl=bar(X,'stacked','BarWidth',.2);

CM=[231,98,84;239,138,71;247,170,88;255,208,111;255,230,183;
    170,220,224;114,188,213;82,143,173;55,103,149;30,70,110]./255;
for i=1:10
    barHdl(i).FaceColor=CM(i,:);
end
% 以下是生成连接的部分
hold on;axis tight 
yEndPoints=reshape([barHdl.YEndPoints]',length(barHdl(1).YData),[])';
zeros(1,length(barHdl(1).YData));
yEndPoints=[zeros(1,length(barHdl(1).YData));yEndPoints];
barWidth=barHdl(1).BarWidth;
for i=1:length(barHdl)
    for j=1:length(barHdl(1).YData)-1
        y1=min(yEndPoints(i,j),yEndPoints(i+1,j));
        y2=max(yEndPoints(i,j),yEndPoints(i+1,j));
        if y1*y2<0
            ty=yEndPoints(find(yEndPoints(i+1,j)*yEndPoints(1:i,j)>=0,1,'last'),j);
            y1=min(ty,yEndPoints(i+1,j));
            y2=max(ty,yEndPoints(i+1,j));
        end
        y3=min(yEndPoints(i,j+1),yEndPoints(i+1,j+1));
        y4=max(yEndPoints(i,j+1),yEndPoints(i+1,j+1));
        if y3*y4<0
            ty=yEndPoints(find(yEndPoints(i+1,j+1)*yEndPoints(1:i,j+1)>=0,1,'last'),j+1);
            y3=min(ty,yEndPoints(i+1,j+1));
            y4=max(ty,yEndPoints(i+1,j+1));
        end
        fill([j+.5.*barWidth,j+1-.5.*barWidth,j+1-.5.*barWidth,j+.5.*barWidth],...
            [y1,y3,y4,y2],barHdl(i).FaceColor,'FaceAlpha',.4,'EdgeColor','none');
    end
end
% 图窗和坐标区域修饰
set(gcf,'Position',[200,100,720,420]);
ax=gca;box off
ax.YLim=[0,1];
ax.XMinorTick='on';
ax.YMinorTick='on';
ax.LineWidth=.8;
ax.FontName='Cambria';
ax.FontSize=11;
ax.TickDir='out';
ax.XTickLabel='10','20','30','40','50','60','70','80';
ax.XLabel.String='Decimals';
ax.YLabel.String='Proportion';
ax.XLabel.FontSize=13;
ax.YLabel.FontSize=13;
ax.Title.String='Area Chart of Proportion — 10-80 digits';
ax.Title.FontSize=14;
% 绘制图例并修饰
lgdHdl=legend(barHdl,num2cell('0123456789'));
lgdHdl.NumColumns=5;
lgdHdl.FontSize=11;
lgdHdl.Location='southeast';

---

5 有向弦图

需要用到这篇的工具函数：https://blog.csdn.net/slandarer/article/details/127992509

统计前后相挨着数字组数量，就比如某两位小数组合14比较多，那么由1指向4的箭头就会更宽一点：

% 构建连接矩阵
dataMat=zeros(10,10);
Pi=getPi(1001);
for i=1:1000
    dataMat(Pi(i)+1,Pi(i+1)+1)=dataMat(Pi(i)+1,Pi(i+1)+1)+1;
end

BCC=biChordChart(dataMat,'Arrow','on','Label',num2cell('0123456789'));
BCC=BCC高分悬赏 限三天 特急matlab  利用凯泽窗函数设计一个带通滤波器

 pi/4QPSK调制解调原理
 pi/4QPSK调制解调原理
 matlab如何标注出点的坐标值
 Matlab中w=(-4:0.001:4)*pi表示啥意思
 MATLAB实现高斯混合分布的EM算法及二维时概率密度曲面置信椭圆绘制