本博客讨论了以下问题:
采样点数N以及采样频率Fs对信号频谱的影响?
最大频谱所对应频率?
产生频率为20Hz和20.02的两个余弦信号的和信号,并分析信号加窗对频率分辨率的影响?
结果1:
结果2:
结果3:
结果4:
结果5:
代码:
clear,clc,close all;
t_start=tic;
%数字信号处理第二次作业
tl=0.5;
f1=20;%频率
f2=20.02;
fs=60;
N=60;
O0=0;
cnt1=1;
fs=50;
for N=10:30:1000
disp(strcat(strcat('第',num2str(cnt1)),'次运算'));
Amp_max1(cnt1,:)=spect_function(f1,fs,N,O0,tl);
cnt1=cnt1+1;
end
close all;
N=200;
cnt2=1
for fs=50:50:2500
disp(strcat(strcat('第',num2str(cnt2)),'次运算'));
Amp_max2(cnt2,:)=spect_function(f1,fs,N,O0,tl);
cnt2=cnt2+1;
end
close all;
cnt1=1;
fs=50;
for N=500:500:5000
disp(strcat(strcat('第',num2str(cnt1)),'次运算'));
Amp_max3(cnt1,:)=Copy_of_spect_function(f1,f2,fs,N,O0,tl);
cnt1=cnt1+1;
end
close all
N=60000;
cnt2=1
for fs=50:50:1000
disp(strcat(strcat('第',num2str(cnt2)),'次运算'));
Amp_max4(cnt2,:)=Copy_of_spect_function(f1,f2,fs,N,O0,tl);
cnt2=cnt2+1;
end
close all;
t_end=toc(t_start);
disp(strcat('程序的运行时间是:',num2str(t_end)));
被调:
function [ Amp_max ] = spect_function( f1,fs,N ,O0,t1)
n=0:N-1;
x=cos(2*pi*f1*n/fs+O0);
figure;
t=0:0.0001:(N-1)/fs;
subplot(311);
plot(t,cos(2*pi*f1*t+O0));title('cos(2*pi*f1*t+O0)');grid on;
subplot(312)
%f = (0:N-1)*fs/N;
% plot(n*fs/N,x,'*');hold on;plot(n*fs/N,x)
plot(n/fs,x,'*');hold on;plot(n/fs,x)
xlabel('t/s');
title(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
grid on;
X = fftshift(fft(x./N));
f = (-N/2:N/2-1)*(fs/N);
subplot(313)
plot(f,abs(X));
hold on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
title(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
grid on;
hold off;
Amp_max=[f(X_f_max_index),X_f(X_f_max_index)];
Amp_max=[f1,fs,N,Amp_max];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=0;
while N1<7
Amp_max=[Amp_max,inf];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=1;
end
if N1>7
Amp_max=Amp_max(1:7);
Amp_max(1,8)=-1;
end
disp(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
disp(strcat(strcat('幅度最大值对应的频率是:',num2str(f(X_f_max_index))),strcat(' 幅度最大值是:',num2str(X_f(X_f_max_index)))));
s1=sum(x.^2);
s2=sum(abs(fft(x)).^2)/N;
Amp_max(1,9)=s1;
Amp_max(1,10)=s2;
flag=0;
if abs(s1-s2)<=10^-5
flag=1;
disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'满足帕塞瓦尔定理'));
else
disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'不满足帕塞瓦尔定理'));
end
Amp_max(1,11)=flag;
disp('----------------------------------');
if(length(X_f_max_index)~=2)
if length(X_f_max_index)>2
X_f_max_index=X_f_max_index(1:2);
end
while length(X_f_max_index)<2
X_f_max_index=[X_f_max_index,inf];
end
end
Amp_max=[Amp_max,X_f_max_index];
pause(t1);
end
被调2:
function [ Amp_max ] = Copy_of_spect_function( f1,f2,fs,N,O0,t1 )
n=0:N-1;
x=cos(2*pi*f1*n/fs+O0)+cos(2*pi*f2*n/fs+O0);
figure;
t=0:0.0001:(N-1)/fs;
subplot(311);
plot(t,cos(2*pi*f1*t+O0)+cos(2*pi*f2*t+O0));title('原始信号cos(2*pi*f1*t+O0)+cos(2*pi*f2*t+O0)');grid on;
subplot(312)
plot(n/fs,x,'*');hold on;plot(n/fs,x)
xlabel('t/s');
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
grid on;
X = fftshift(fft(x./N));
f = (-N/2:N/2-1)*(fs/N);
subplot(313)
plot(f,abs(X));
hold on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
grid on;
hold off;
figure;
subplot(311)
plot(f,abs(X));
hold on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
grid on;
hold off;
subplot(312)
plot(f,abs(X));
hold on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
axis([-20.1 ,-19.9,-inf inf]);
title('负频谱');
grid on;
hold off;
subplot(313)
plot(f,abs(X));
hold on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
axis([19.9 ,20.1,-inf inf]);
title('正频谱');
grid on;
hold off;
Amp_max=[f(X_f_max_index),X_f(X_f_max_index)];
Amp_max=[f1,fs,N,Amp_max];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=0;
while N1<7
Amp_max=[Amp_max,inf];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=1;
end
if N1>7
Amp_max=Amp_max(1:7);
Amp_max(1,8)=-1;
end
disp((strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))))));
disp(strcat(strcat('幅度最大值对应的频率是:',num2str(f(X_f_max_index))),strcat(' 幅度最大值是:',num2str(X_f(X_f_max_index)))));
s1=sum(x.^2);
s2=sum(abs(fft(x)).^2)/N;
Amp_max(1,9)=s1;
Amp_max(1,10)=s2;
flag=0;
if abs(s1-s2)<=10^-5
flag=1;
disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'满足帕塞瓦尔定理'));
else
disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'不满足帕塞瓦尔定理'));
end
Amp_max(1,11)=flag;
disp('----------------------------------');
pause(t1);
end