待修改7777

clc
clear 
%批量读取一个文件夹txt数据
cell_type = 'pc';
namelist=dir('D:\zhuomian\shuju\pc\*.txt'); 
l=length(namelist);
for i=1:l
    file_name=strcat(['D:\zhuomian\shuju\pc\',namelist(i).name]);
    Data{i}=textread(file_name); %不同的文件此处对应的函数不同，如excel文件则是xlxread(filename)；
end


% 获取第一个文件的大小来定义 data
firstFile = textread(strcat(['D:\zhuomian\shuju\pc\', namelist(1).name]));
dataSize = size(firstFile);
data1 = zeros(l, dataSize(1));
% 初始化 data1
for i = 1:l
 data = Data{i};
data(:,2) = smoothdata(data(:,2),'sgolay');
 signal = data(:,2);
 data(:,i) = wdenoise(signal,5,'DenoisingMethod','BlockJS'); %%采用小波去噪函数不在单独平滑 
[xd1,cxd1,lxd1] = wden(signal,'minimaxi','s','sln',3,'db3');
 data(:,2) = xd1;
data(:,2) = smoothdata(data(:,2),'gaussian',5);%平滑
data(:,1) = round(data(:,1));%Raman shift 取整
 source_data = data(:,2);
MappedData = mapminmax(source_data, -1, 1);% 最值归一化
normalized_data = zscore(source_data); %标准化
 data1(i,:) = normalized_data;
end 
figure
plot(data(:,1),mean(data1(1:10,:)+1),'g',LineWidth=3);hold on
plot(data(:,1),mean(data1(11:20,:)+1),'r',LineWidth=3);hold on
plot(data(:,1),mean(data1(21:30,:)+1),'b',LineWidth=3);hold on
xlabel('Raman shift/cm^-1',LineWidth=3)
 ylabel('Intenstiy',LineWidth=3)
title('Raman spectroscopy',LineWidth=3)
legend('PANC-NC','PANC-OE','PANC-1')

% % 修改图例，将图例的标签设置为文件名。
% legend(namelist(1:l).name);
% 设置是否输出图像的标志
outputImage = false; % 设置为 true 以输出图像

if outputImage
    % 修改图例，将图例的标签设置为文件名
    legend(namelist(1:l).name);

% 绘制数据
figure;

for i = 1:l
    plot(Data{i}(:,1), data1(i,:), 'LineWidth', 1.5);
    hold on;
end

end


% 1. 峰对齐
alignedData = zeros(size(data1));
referenceSpectrum = data1(1, :); % 将第一个光谱作为参考

for i = 1:l
    % 使用交叉相关进行峰对齐
    [xc, lag] = xcorr(referenceSpectrum, data1(i, :));
    [~, idx] = max(xc);
    alignmentShift = lag(idx);

    % 对光谱进行平移以进行对齐
    alignedData(i, :) = circshift(data1(i, :), [0, alignmentShift]);
end

% 2. 特征提取（示例：提取每个光谱的均值和标准差）
meanValues = mean(alignedData, 2); % Mean value for each spectrum
stdValues = std(alignedData, 0, 2); % Standard deviation for each spectrum

% 3. 特征峰识别（示例：使用findpeaks函数进行简单的峰检测）
peakIndices = cell(l, 1);
peakValues = cell(l, 1);

for i = 1:l
    [peakValues{i}, peakIndices{i}] = findpeaks(alignedData(i, :), 'MinPeakHeight', 0.5, 'MinPeakDistance', 10);
end

% 4. 保存每个光谱的峰值位置和强度
peakPositions = zeros(l, 1);
peakIntensities = zeros(l, 1);

for i = 1:l
    % 使用 findpeaks 找到峰值和对应的位置、强度
    [peaks, locs] = findpeaks(alignedData(i, :));

    % 如果有峰值，选择位置和强度最大的峰
    if ~isempty(peaks)
        [~, index] = max(peaks);
        peakPositions(i) = locs(index);
        peakIntensities(i) = peaks(index);
    else
        peakPositions(i) = NaN;  % 如果没有峰值，标记为 NaN
        peakIntensities(i) = NaN;
    end
end

% 绘制直方图
figure;
histogram(peakPositions, 'BinWidth', 1);  % 可以根据需要调整 BinWidth
title('峰值位置直方图');
xlabel('峰值位置');
ylabel('频数');

% 统计分析
meanDifference = mean(peakPositions);
medianDifference = median(peakPositions);
stdDevDifference = std(peakPositions);

fprintf('平均峰值位置差异: %.2f\n', meanDifference);
fprintf('中位数峰值位置差异: %.2f\n', medianDifference);
fprintf('峰值位置差异标准差: %.2f\n', stdDevDifference);

% 5. 特征区域选择
selectedRegionIndices = find(data1(1, :) > 0.5); % 选择数据大于某个阈值的区域

% 6. 峰强度、位置和宽度提取
peakIntensity = cell(l, 1);
peakPosition = cell(l, 1);
peakWidth = cell(l, 1);

for i = 1:l
    % 使用 findpeaks 函数找到峰值的强度、位置和宽度
    [peaks, locs, widths] = findpeaks(alignedData(i, :), 'MinPeakHeight', 0.5, 'MinPeakDistance', 10);

    % 保存峰值的强度、位置和宽度
    peakIntensity{i} = peaks;
    peakPosition{i} = locs;
    peakWidth{i} = widths;
end


% 7. 可视化结果
figure;

% 绘制峰对齐后的数据
for i = 1:l
    plot(Data{i}(:,1), alignedData(i, :), 'LineWidth', 1.5);
    hold on;
end

% 在图上标记特征峰
for i = 1:l
    plot(Data{i}(peakIndices{i}, 1), alignedData(i, peakIndices{i}), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
end

legend(namelist(1:l).name);


% 显示其他结果，根据需要自行调整
disp('均值:');
disp(meanValues);

disp('标准差:');
disp(stdValues);

disp('峰值位置:');
disp(peakPosition);

disp('峰值强度:');
disp(peakIntensity);

disp('峰值宽度:');
disp(peakWidth);

待修改