基于MATLAB光谱数据预处理.rar

clear all; clc; load('A.mat'); load('RAW.mat'); %400*255 255列是样本值 400为样本数 x=A(1,:); %x为张柳师兄，前300个为好种子，最后100个为得了病的种子，赤霉病是白的，所以反射率会变大 y=RAW(:,1:254); %y为自己提取的感兴趣区域，但只要第一行波长的信息 z=RAW(:,255); %b1=A(:,:); %y1=b1(2:6,:); %[m n]=size(y);返回矩阵y的尺寸信息, 并存储在m、n中。其中m中存储的是行数,n中存储的是列数。 %% 两种种子平均光谱对比 % haode=y(1:300,:); % huaide=y(301:400,:); % haode=mean(haode); % huaide=mean(huaide); % figure(1); % plot(x,haode); % axis([866.4,1701,0,1.0]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('两种种子对比'); % hold on; % plot(x,huaide); % legend('好的种子','坏的种子'); %% 画原始光谱图和平均光谱图 % figure(2); % plot(x,y); % axis([866.4,1701,0,1.0]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('张柳师兄数据'); %% 画四组数据的平均波形图 % y1=mean(y(1:100,:)); % figure(3); % plot(x,y1); % axis([866.4,1701,0,1.0]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('张柳师兄数据四组平均光谱'); % hold on; % y2=mean(y(101:200,:)); % plot(x,y2); % y3=mean(y(210:300,:)); % plot(x,y3); % y4=mean(y(301:400,:)); % plot(x,y4); %% 验证平滑的作用 % y=Y+(0.4-rand(400,254));%Y就是400*254的反射率 % [smdata] = SMOOTH(y,5); % [smdata1] = SMOOTH(y,11); % [smdata2] = SMOOTH(y,21); % figure(4); % plot(x,Y(1,:),'k-',x,y(1,:),'m-x',x,smdata(1,:),'c-o',x,smdata1(1,:),'b-p',x,smdata2(1,:),'r-+'); % legend('原始值','加入随机噪声','5点平滑','11点平滑','21点平滑'); % axis([866.4,1701,-0.5,1.5]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); %% 光谱数据预处理，先进行降噪，所谓降噪就是求导，卷积等方法来滤波，所以降噪完的图像和原图像有一定差距，但是数据的趋势等不发生变化 %先用SNV(标准正态变量变换）降噪，看一下降噪完的光谱图 % y_SNV=snv(y(:,1:254)); % figure(4); % plot(x,y_SNV); % axis([866.4,1701,-2,2.0]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('SNV预处理后的光谱图像'); %用SNV降噪完后再用SG卷积平滑来继续降噪 % y_SNV_SG=savgol(y_SNV,15,3,1); % figure(5); % plot(x,y_SNV_SG); % axis([866.4,1701,-0.2,0.2]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('SNV之后在进行卷积平滑预处理后的光谱图像'); %% 用MSC（多元散射校正）降噪，选的最优波形图是所有数据的平均，看一下降噪完的光谱图 %y_MSC=msc(y,y_m); % figure(6); % plot(x,y_MSC); % axis([866.4,1701,0.3,0.9]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('MSC预处理后的光谱图像(外部参照是所有光谱的平均值）'); %上面的外部参照选的是所有光谱的平均值，下面的外部参照选前三百个好的种子的平均值 % y_m_1=mean(y(1:300,:)); % y_MSC_1=msc(y,y_m_1); % figure(7); % plot(x,y_MSC_1); % axis([866.4,1701,0.2,0.9]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('MSC预处理后的光谱图像(外部参照是前三百个光谱的平均值）'); %% 归一化Normalize来进行预处理 % y_nor=normaliz(y); % figure(8); % plot(x,y_nor); % axis([866.4,1701,0,0.1]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('归一化预处理后的光谱图像'); %% SG卷积平滑法及求导 % y_SG=savgol(y,3);%窗口大小必须为奇数，它的原理是最小二乘法来进行平滑，所有需要奇数个点来进行平滑后替代中间多出来的点 % %但是这步本身没有意义，y_SG和原矩阵一模一样 % y_SG_1=savgol(y_SG,5,3,1);%光谱数据，平滑点数，多项式，一阶求导 % figure(9); % plot(x,y_SG_1); % axis([866.4,1701,-0.05,0.05]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('SG卷积平滑法预处理后的光谱图像'); %% 小波变换预处理光谱 % y_DWT=DWT(y); % figure(10); % plot(x,y_DWT); % axis([866.4,1701,0,1]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('小波变换预处理后的光谱图像'); %% 微分预处理光谱，一般不用 % y1st=diff(y,1);%X?为输入光谱矩阵，X1st?为输出一阶导数光谱矩阵 % figure(11); % plot(x,y1st); % axis([866.4,1701,-1,1]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('一阶微分预处理后的光谱图像'); %y2st=diff(y,2);%X?为输入光谱矩阵，X2st?为输出二阶导数光谱矩阵 %% 窗口平滑处理（smooth） % [smdata] = SMOOTH(y,3);%或者[mdata]=nirmaf(y,3)，窗口必须是奇数,窗口越大，失真越严重 % figure(20); % plot(x,smdata); % axis([866.4,1701,0,1]); % xlabel('Wavelength/nm','FontSize',12); % ylabel('Refelctance/%','FontSize',12); % title('平滑预处理后的光谱图像'); %% 直接正交信号校正（DOSC） [Z,W,P,T] = dosc(y,z,2,1E-3); figure(21); subplot(121); plot(x,y); title('原始数据') subplot(122); plot(x,Z); title('直接正交校正后的数据') Xcorr = corrcoef([y z]); Zcorr = corrcoef([Z z]); figure(22); plot(Xcorr(1:end-1,end),'b-'); hold on plot(Zcorr(1:end-1,end),'r:') %% 连续投影算法（SPA）提取特征波长，Xcal是训练集，无标签；Xval是验证集，无标签；ycal是训练集的标签；yval是验证集的标签。先把标签给放回来 % z=RAW(:,255); % y_SPA_RAW=[y_MSC,z]; % % % %第一种标签，每100行一个标签 % [m n] = size(y_SPA_RAW); %m样本数，n维数 % R = randperm(m); %1到m这些数随机打乱得到的一个随机数字序列作为索引 % RAWtest = y_SPA_RAW(R(1:100),:); %以索引的前100行作为测试样本RAWtest % R(1:100) = []; %空中括号表示删去 % RAWtraining = y_SPA_RAW(R,:); %剩下的数据作为训练样本RAWtraining % RAW_c=RAWtraining(:,1:n-1); % RAW_v=RAWtest(:,1:n-1); % RAWc=RAWtraining(:,n); % RAWv=RAWtest(:,n); % Xcal=RAW_c; % ycal=RAWc; % Xval=RAW_v; % yval=RAWv; % % [m n] = size(y_SPA_RAW); %m样本数，n维数 % R = randperm(m); %1到m这些数随机打乱得到的一个随机数字序列作为索引 % RAWtest = y_SPA_RAW(R(1:100),:); %以索引的前100行作为测试样本RAWtest % R(1:100) = []; %空中括号表示删去 % RAWtraining = y_SPA_RAW(R,:); %剩下的数据作为训练样本RAWtraining % [m1 n1]=size(RAWtraining); % R1 = randperm(m1); % RAWtest1 = RAWtraining(R1(1:75),:); % R1(1:75) = []; % RAWtraining1 = RAWtraining(R1,:); % % RAW_c=RAWtraining1(:,1:n1-1); % RAW_v=RAWtest1(:,1:n1-1); % RAWc=RAWtraining1(:,n1); % RAWv=RAWtest1(:,n1); % Xcal=RAW_c; % ycal=RAWc; % Xval=RAW_v; % yval=RAWv; %% 第二种标签，前300行为1，后100行为2 % % load('RAW_another.mat'); % % [m n] = size(RAW_another); %m样本数，n维数 % % R = randperm(m); %1到m这些数随机打乱得到的一个随机数字序列作为索引 % % RAW_anothertest = RAW_another(R(1:100),:); %以索引的前100行作为测试样本RAW_anothertest % % R(1:100) = []; %空中括号表示删去 % % RAW_anothertraining = RAW_another(R,:); %剩下的数据作为训练样本RAW_anothertraining % % RAW_another_c=RAW_anothertraining(:,1:n-1); % % RAW_another_v=RAW_anothertest(:,1:n-1); % % RAW_anotherc=RAW_anothertraining(:,n); % % RAW_anotherv=RAW_anothertest(:,n); % % Xcal=RAW_another_c; % % ycal=RAW_anotherc; % % Xval=RAW_another_v; % % yval=RAW_anotherv; % % % % % % % RAW_spa=spa(Xcal,ycal,Xval,yval,5,20,1);%可以调用spa函数，特征波长选择最少5个，最多不限，选择自动校准,自动校准（autoscaling）选择1，否则选择0 % % m_min=5; % m_max=50; % % N = size(Xcal,1); %Xcal的行数，也就是样本数 % K = size(Xcal,2); %Xcal的列数，也就是维数 % % % %% Phase 1: Projection operations for the selection of candidate subsets选择候选子集的投影运算 % % normalization_factor = std(Xcal);%求Xcal的标准偏差，也就是标准化因子，化数据为一行。标准化的作用：标准化之后使不同的特征具有相同的尺度 % % for k = 1:K % x = Xcal(:,k); % Xcaln(:,k) = (x - mean(x)) / normalization_factor(k); %每一列数据减去平均值除以标准差，就是标准差标准化（