傅里叶频谱分析

做音频的谱图分析，工具为MATLAB，不多说，直接上代码

clear all;
close all;
clc;
diroutput=dir('E:\000\0deep\AccelScansComponents\train\*.txt')

index_cell_train={'G1SquaredAluminumMesh';'G2GraniteTypeVeneziano';'G3AluminumPlate';'G4Bamboo';'G5SolidRubberPlateVersion1';'G6Carpet';'G7FineFoamVersion1';'G8Cardboard';'G9Jeans'};

in=[1;2;3;4;5;6;7;8;9;10];
Fs=10000;
% load mycolor
Fs=10000;
window=hamming(512);
noverlap=128;
nfft=512;
bi=1;
%fid=fopen('train.txt','a');

for j=1:9

for i=1:9
    p=in(i);
bi
%ax=textread('E:\000\0deep\AccelScansComponents\train\G1SquaredAluminumMesh_Movement_X_train1.txt','%n','whitespace','');
ax=textread(['E:\000\0deep\AccelScansComponents\train\',index_cell_train{i},'_Movement_X_train',num2str(p)','.txt'],'%n','whitespace','');
ay=textread(['E:\000\0deep\AccelScansComponents\train\',index_cell_train{i},'_Movement_Y_train',num2str(p)','.txt'],'%n','whitespace','');
az=textread(['E:\000\0deep\AccelScansComponents\train\',index_cell_train{i},'_Movement_Z_train',num2str(p)','.txt'],'%n','whitespace','');


ax1= fft(ax);
ay1= fft(ay);
az1= fft(az);
as=sqrt(((abs(ax1)).^2)+((abs(ay1)).^2)+((abs(az1)).^2));% frequency-domain magnitude of our new DFT321 signal
Im=imag(ax1+ay1+az1);
Re=real(ax1+ay1+az1);
B=Im./Re;
AJ=atan(B);
AS=complex(as.*cos(AJ),as.*sin(AJ));
IAS=ifft(AS');
a=lpc(IAS,400);
est_IAS=filter([0 -a(2:end)],1,IAS);
Sample_Movement_Acc_Train{bi,:}=IAS;

hCircle=figure('Visible','off');%guanbichuangkou
hold on;
[B, F, T, P] = spectrogram(est_IAS,window,noverlap,0:1:512,Fs,'yaxis');%spectrogram(X,256,250,256,2000)yaxis是切换时间对频率;
%mymode=gray;
  %mymode=mymode(64:-1:1,:);
 % colormap(mymode);
%B = mapminmax(B,0,1)
% imshow(P,[])
F = F(1:256)
T = T(1:43)
  Z = imagesc(T,F,log10(abs(B)));
  image = imresize(Z,[128 128])
  %imshow(image)
  %imwrite(image,'E:/000/0deep/fit2/i.jpg')
  
  set(gca,'YDir','normal');%a4=colormap(gray);%让Y轴顺序正常
axis on%time an frequence
saveas(hCircle,['E:/000/0deep/fit2/',index_cell_train{j},'train',num2str(p),'.jpg'],'jpg');
%fprintf(fid,'%s %d\n',[index_cell_train{j},'train',num2str(p),'.jpg'],9);

 bi=bi+1;
end

end

fs = 10000;
%t = 1.65
t = 4

diroutput=dir('E:\000\0deep\AccelScansComponents\train\*.txt')

index_cell_train={'G1SquaredAluminumMesh';'G2GraniteTypeVeneziano';'G3AluminumPlate';'G4Bamboo';'G5SolidRubberPlateVersion1';'G6Carpet';'G7FineFoamVersion1';'G8Cardboard';'G9Jeans'};
in=[1;2;3;4;5;6;7;8;9;10];
bi = 1
% for j=1:9
for pp=1:9
    p=in(pp)
bi
%ax=textread('E:\000\00zj\AccelScansComponents\train\G1SquaredAluminumMesh_Movement_X_train1.txt','%n','whitespace','');
ax=textread(['E:\000\00zj\accels\',index_cell_train{pp},'_Movement_X_train',num2str(p)','.txt'],'%n','whitespace','');
ay=textread(['E:\000\00zj\accels\',index_cell_train{pp},'_Movement_Y_train',num2str(p)','.txt'],'%n','whitespace','');
az=textread(['E:\000\00zj\accels\',index_cell_train{pp},'_Movement_Z_train',num2str(p)','.txt'],'%n','whitespace','');

[b,a]=butter(5,10/(10000/2),'high');
ax=filter(b,a,ax);
ay=filter(b,a,ay);
az=filter(b,a,az);

ax1= fft(ax);
ay1= fft(ay);
az1= fft(az);
as=sqrt(((abs(ax1)).^2)+((abs(ay1)).^2)+((abs(az1)).^2));% frequency-domain magnitude of our new DFT321 signal

Im=imag(ax1+ay1+az1);
Re=real(ax1+ay1+az1);
B=Im./Re;
AJ=atan(B);
AS=complex(as.*cos(AJ),as.*sin(AJ));
IAS=ifft(AS);
a=lpc(IAS,400);
est_IAS=filter([0 -a(2:end)],1,IAS);
y=est_IAS
%y = sin(128*pi*t) + sin(256*pi*t);      

figure;
win_sz = 512;
han_win = hanning(win_sz);      % 选择海明窗

nfft = win_sz;
nooverlap = 128;
[S, F, T] = spectrogram(y, window, nooverlap, nfft, fs);
% T = T(1:3)

hCircle=figure('Visible','on');%guanbichuangkou
imagesc(T, F, log10(abs(S)))
set(gca, 'YDir', 'normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('short time fourier transform spectrum')
saveas(hCircle,['E:/000/0deep/fit2/',index_cell_train{pp},'train',num2str(p),'.jpg'],'jpg');

bi=bi+1;
end
% end