写在开头：临近期末，正在复习Matlab课程，顺便将自己的感悟写下，仅作记录，如有错误欢迎指正。

III.Siganlanalyse

1.对坐标数据进行处理

目标：给一个包含(x，y，z)信息的txt文件，从Matlab中读取，并将它plot出来。

S=importdata('文件名')H=S.data.*1e-3; %s里的数据单位保持一致(毫米转化成米)L=4e-3 %数据的长度 从txt文件中得知，或者使用nomal()X=size(H,1)Y=size(H,2)%Size函数获得的是目标的维度信息，如[3,4]X=linspace(0,L,X);Y=linspace(0,L,Y);%Y从1->L 每隔Y个间距，一个点surf(X,Y,H); %将 三维信息 展示出来zlim(min(min(H)), max(max(H))); %H是峰值，所以要对z轴的高度进行限制，H是一个矩阵max(H)是矩阵中所有列最大项，max(max(H))为所有列最大项中的最大项，即为H矩阵中最大项。title('Topographie');%添加标题xlabel('x-position[m]');ylabel('y-position[m]');zlabel('Hohe[m]'); %对坐标轴添加信息

Ubung练习1

% Programmbeschreibung: Ubung zur Vorlesung 9% Autor:Jiangshunxin% Datum:24.1.2021clear allclose allclc%Pragramm Anfang%%S=importdata("data_time.txt");t=S.data(:,1);L=S.data(:,2:5);figureplot(t,L);title('U9');xlabel('x-position[t]');ylabel('y-position Mess');legend({'L1','L2','L3','L4'},'location','best')xlim([min(t) max(t)]);grid on;box on;save('data_time','t','L');saveas(gcf,'U9A1.fig');saveas(gcf,'U9A1.png');%%%Programm Ende

t=s.data()-读取数据,

gcf-获取当前图窗句柄

xlim[min() max()]限制x轴显示区域

box on 显示坐标区框轮廓

U9数据处理图像

2.对信号数据进行处理

1.信号的直接处理

• 使用离散傅里叶变换

[y,Fs]=audioread('文件名'); % y是采样数据(时间数据向量)；Fs是采样率N=momel(y); %获取采样数据的采样点个数Yc=fft(y)/N; %计算FFT，fft是N个测量点的信息，Yc是包含幅值和相位移的复数信息的矩阵Y=2*abs(Yc(1:N/2+1)); %Yc(1:N/2+1)去除Yc信息矩阵中复数部分，abs()将Yc中负数部分去除，故需要*2Y(1)=Y(1)/2; %由于上下对称性，确定幅值/2fB=Fs/2;%信号的带宽，DFT中只有实数值，由于共轭所以/2f=fB*linspace(0,1,N/2+1); %df=fs/N =fB/N/2, 此处f为从fB->0间距递减的值，为x轴坐标。plot(f,Y);title('Frequenzspektrum des Quietschens');ylabel('Amplitude[-]');xlabel('Frequenz[Hz]');

• 带宽fB是采样率(Fs)的一半

• 带宽Fb包含(Nf=N/2+1)个频率线，它们将Fb均分。

2.窗函数解决频谱泄露

函数：使用汉明窗

win=windows(@handle,N)->win = hann(N)

• win:权重因子

• N:采样点数

使用匿名函数

hannwin=@(N)0.5*(1+cos((2*pi*linspace(-N/2,N/2-1,N))/N))；function win = hannwin(N)nWin = linspace(-N/2,N/2-1,N);win=0.5*(1+cos((2*pi*nWin)/N));end 

优化后的信号处理代码

[y,Fs]=audioread('文件名'); % y是采样数据(时间数据向量)；Fs是采样率N=numel(y); %获取采样数据的采样点个数win=hann(N);y=y.*win;Yc=fft(y)/N; %计算FFT，fft是N个测量点的信息，Yc是包含幅值和相位移的复数信息的矩阵Y=2*abs(Yc(1:N/2+1)); %Yc(1:N/2+1)去除Yc信息矩阵中复数部分，abs()将Yc中负数部分去除，故需要*2Y(1)=Y(1)/2; %由于上下对称性，确定幅值/2fB=Fs/2;%信号的带宽，DFT中只有实数值，由于共轭所以/2f=fB*linspace(0,1,N/2+1); %df=fs/N =fB/N/2, 此处f为从fB->0间距递减的值，为x轴坐标。hannwin=@(N)0.5*(1+cos((2*pi*linspace(-N/2,N/2-1,N))/N))；function win = hannwin(N)nWin = linspace(-N/2,N/2-1,N);win=0.5*(1+cos((2*pi*nWin)/N));endplot(f,Y);title('Frequenzspektrum des Quietschens');ylabel('Amplitude[-]');xlabel('Frequenz[Hz

3.光谱图

• 使用快速傅里叶变换对无穷周期信号定义

公式：N=[1+(ns-1)(1-overlap)]*L

N:信号的整体长度

L:单个片段的长度

ns：片段数

overlap：单个片段的重叠百分比

nsec=500; %片段数overlap = 0.5; %定义重叠百分比L=floor(N/(1+(nsec-1)*(1-overlap))); %floor将每个元素四舍五入得到整数L=L-mod(L,2);  %精确计算L的长度start = floor(linspace(1,N-L,nsec));ende = start+L-1; Amp = zeros(L/2+1,nsec); %准备好存储的0矩阵win = hann(L); %汉明窗函数for idx = 1:nsec % 使用FFT计算  y_idx = y(start(idx):ende(idx)); % 选择段落  y_idx = y_idx.*win*L/sum(win); % 窗函数并修正峰值  Yc = fft(y_idx)/L; % 使用快速傅里叶变换计算  Amp(:,idx) = 2*abs(Yc(1:L/2+1));   Amp(1,idx) = Amp(1,idx)/2; %计算峰值endfigure imagesc(t,f/1000,Amp); % 补充颜色%imagesc(t,f/1000,log10(Amp); %使用Log函数让图谱颜色区更深c = colorbar; % Drehen der Ansicht view(0,-90); % 旋转视图title('Spektrogramm "Quietschen.wav"'); xlabel('Zeit in s'); ylabel('Frequenz in kHz'); ylabel(c,'Amplitude in -'); %caxis([0 0.01]); %添加刻度%%%%%%%%%%%%%%%%%%%%%figure spectrogram(y,hann(L),L/2,L,Fs, 'yaxis'); title('Spektrogramm mit Toolbox'); %使用Toolbox函数直接输出图像，y-数据向量，win-权重，L/2-重叠数据点的数量，%L-快速傅里叶变换中数据点的数量，Fs-采样率，'yaxis'-频率在y轴显示

4.数字滤波

1. 计算过滤因子(Butterworth Filter)

[b,a] = butter(ord, Wn, 'ftype');

[b,a]-Filterkoeffizienten,ord-过滤阶数 ,Wn-Normierte Grenzfrequenz,'ftype'-过滤类型

Normierte Grenzfrequenz:

Wn=fg/fb=fg/(0.5*Fs)

• fg:Grenzfrequenz,限制频率

• fb:Signalbandbreite,信号带宽

• Fs:Abtastfrequenz,采样率

Filtertypen:

'high'->Hochpass-高通，'low'->低通，

'stop'->Bandsperre Wn=[w1,w2], '-'->Bandpass Wn=[w1,w2]

1. 计算 ord和Wn 阶数和过滤因子

[ord,Wn] = buttord(Wp,Ws,Rp,Rs);

• Rs: Min.Sperrdämpfung in dB/最小抑制波段衰减

• Rp: Durchgangsdämpfung in dB/通过波段

• Ws: Normierte Sperrfrequenz/标定抑制频率

• Wp: Normierte Durchlassfrequenz/标定通过频率

1. 描绘分析过滤

freqz(b,a,n,Fa);

b,a->Filters 的因子，n-Frequenzauflösung频率选择，Fa-采样比率-类似50%，

[h,tf]=freqz(b,a,n,Fa);

输出峰值和相位曲线在复数矩阵 h和对应频率 tf

freqz(b,a,1000,500)，n=1000,Fa=500

1. 使用Matlab中的过滤函数

• fftfilt()->FFT-based FIR filtering

• filtfilt()->Zero-phase digital filtering

...[ord,Wn] = buttord(5000/0.5/Fs,6500/0.5/Fs,3,40); [b,a] = butter(ord,Wn, 'low'); y = filtfilt(b,a,y); %使用滤波函数...

Ubung 练习2

%Programmbeschreibung:%Ubung9 A2 zur Vorlesung9%Autor:Jiangshunxin%Datum:24.1.2021clear allclose allclc%%%Programm Anfang[y,Fs]=audioread('Sirene.mp3');y=y(:,1); %选择声道N=numel(y);N=N-mod(N,2); %剔除被2除后余数部分，因为要处理N/2y=y(1:N);win = hann(N);yf=y.*win;yf=yf.*N./sum(win);%使用汉明窗修正Yc = fft(y)/N;Ycf = fft(yf)/N;%使用快速傅里叶变换计算频率图谱Amp = 2*abs(Yc(1:N/2+1));Amp(1)=Amp(1)/2;Ampf = 2*abs(Ycf(1:N/2+1));Ampf(1)=Ampf(1);%计算峰值图谱fb=Fs/2;f = fb*linspace(0,1,N/2+1);%对频率分段figureplot(f,Amp,'b-',f,Ampf,'r:');xlabel('Frequenz in Hz');ylabel('Amplitude in -');title('Frequenzspektrum');legend({'ungefensterte Daten','Gefensterte Daten'});grid on; box on;xlim([0 8000]);%调整显示区域从0到8000saveas(gcf,'U9A2.png');saveas(gcf,'U9A2.fig');%%%Programm Ende

U2练习频谱图