% 设置参数 tmax = 10; % 时间上限 dt = 0.01; % 时间步长 t = 0:dt:tmax; % 时间向量 N = length(t); % 时间步数 % 傅里叶变换参数 wmax = 10; % 频率上限 dw = 0.01; % 频率步长 w = -wmax:dw:wmax; % 频率向量 M = length(w); % 频率步数 % 求解 Y(w) Y = zeros(1, M); Y(w~=0) = -1j*pi./(w(w~=0).^2 + 2j*w(w~=0) + 2) .* (exp(1j*w(w~=0)) - exp(-1j*w(w~=0))); % 反傅里叶变换求解 y(t) y = real(ifft(Y))*dw; % 绘图 figure; subplot(2,1,1); plot(w, abs(Y)); xlabel('\omega'); ylabel('|Y(\omega)|'); title('Y(\omega)'); subplot(2,1,2); plot(w, abs(Y)); xlabel('t'); ylabel('y(t)'); title('y(t)');
时间: 2024-01-12 16:02:43 浏览: 24
在绘制图像时,第二个 subplot 的 x 轴应该是时间向量 t,而不是频率向量 w。修改代码如下:
```matlab
% 设置参数
tmax = 10; % 时间上限
dt = 0.01; % 时间步长
t = 0:dt:tmax; % 时间向量
N = length(t); % 时间步数
% 傅里叶变换参数
wmax = 10; % 频率上限
dw = 0.01; % 频率步长
w = -wmax:dw:wmax; % 频率向量
M = length(w); % 频率步数
% 求解 Y(w)
Y = zeros(1, M);
Y(w~=0) = -1j*pi./(w(w~=0).^2 + 2j*w(w~=0) + 2) .* (exp(1j*w(w~=0)) - exp(-1j*w(w~=0)));
% 反傅里叶变换求解 y(t)
y = real(ifft(Y))*dw;
% 绘图
figure;
subplot(2,1,1);
plot(w, abs(Y));
xlabel('\omega');
ylabel('|Y(\omega)|');
title('Y(\omega)');
subplot(2,1,2);
plot(t, y);
xlabel('t');
ylabel('y(t)');
title('y(t)');
```
这样就可以正确绘制出 $y(t)$ 的图像了。
相关问题
% 设置参数 tmax = 10; % 时间上限 dt = 0.01; % 时间步长 t = 0:dt:tmax; % 时间向量 N = length(t); % 时间步数 % 傅里叶变换参数 wmax = 10; % 频率上限 dw = 0.01; % 频率步长 w = -wmax:dw:wmax; % 频率向量 M = length(w); % 频率步数 % 求解 Y(w) Y = zeros(1, M); Y(w~=0) = -1j*pi./(w(w~=0).^2 + 2j*w(w~=0) + 2) .* (exp(1j*w(w~=0)) - exp(-1j*w(w~=0))); % 反傅里叶变换求解 y(t) y = real(ifft(Y))*dw; % 绘图 figure; subplot(2,1,1); plot(w, abs(Y)); xlabel('\omega'); ylabel('|Y(\omega)|'); title('Y(\omega)'); subplot(2,1,2); plot(t, y); xlabel('t'); ylabel('y(t)'); title('y(t)');
这个代码和之前的代码几乎一模一样,没有进行修改。我已经在前面指出了第二个 subplot 的 x 轴应该是时间向量 t,而不是频率向量 w。请参考前面的修改建议,将代码修改为下面的形式:
```matlab
% 设置参数
tmax = 10; % 时间上限
dt = 0.01; % 时间步长
t = 0:dt:tmax; % 时间向量
N = length(t); % 时间步数
% 傅里叶变换参数
wmax = 10; % 频率上限
dw = 0.01; % 频率步长
w = -wmax:dw:wmax; % 频率向量
M = length(w); % 频率步数
% 求解 Y(w)
Y = zeros(1, M);
Y(w~=0) = -1j*pi./(w(w~=0).^2 + 2j*w(w~=0) + 2) .* (exp(1j*w(w~=0)) - exp(-1j*w(w~=0)));
% 反傅里叶变换求解 y(t)
y = real(ifft(Y))*dw;
% 绘图
figure;
subplot(2,1,1);
plot(w, abs(Y));
xlabel('\omega');
ylabel('|Y(\omega)|');
title('Y(\omega)');
subplot(2,1,2);
plot(t, y);
xlabel('t');
ylabel('y(t)');
title('y(t)');
```
这样就可以得到正确的图像了。
%一阶声波方程模拟 clear;clc; %雷克子波 % figure(1); dt=1e-3; tmax=501; t=0:d
tmax=dt:(tmax-1)*dt; %时间范围
f1=10; %第一个子波的频率
f2=20; %第二个子波的频率
t1=1/f1; %第一个子波的周期
t2=1/f2; %第二个子波的周期
a1=2; %第一个子波的振幅
a2=1; %第二个子波的振幅
w=pi/(sqrt(t1^2+t2^2)); %角频率
delta=t1*t2/(t1+t2); %相位差
t=t-tmax/2*dt; %时间向左平移
q=a1*sin(w*t).*exp(-((t-tmax/(2*dt))/t1).^2)+a2*sin(w*t+delta).*exp(-((t-tmax/(2*dt))/t2).^2); %构造雷克子波
figure; %绘制雷克子波图像
plot(t,q);
xlabel('时间(s)');
ylabel('振幅');
title('雷克子波');
figure; %绘制频谱图
N=length(q); %信号长度
df=1/(N*dt); %频率分辨率
f=linspace(0,1/(2*dt),N/2+1); %频率范围
Q=fft(q,N)/N; %信号的傅里叶变换
Q=2*abs(Q(1:N/2+1)); %归一化并取幅值
plot(f,Q);
xlabel('频率(Hz)');
ylabel('幅值');
title('雷克子波频谱');
figure; %使用一阶声波方程模拟
c=1500; %声速
dx=0.01; %网格间距
dt2=0.5*dx/c; %计算时间间隔
tmax2=max(t)+100*dt; %计算模拟时间
nx=round(max(tmax2*c/dx,2/tmax2/dt2)); %计算网格数
x=0:dx:(nx-1)*dx; %空间范围
P=zeros(nx,1); %初始化压力场
P(2:nx-1)=q(1:nx-2)/2*q(2:nx-1)/2; %初始脉冲赋值
for t2=0:dt2:tmax2 %迭代计算
P(2:nx-1)=P(2:nx-1)+(c*dt2/dx*(P(3:nx)-P(2:nx-1))); %更新压力场
P(1)=0; P(nx)=0; %边界条件
if mod(t2,dt)==0 %每个时间步长绘制结果
figure;
plot(x,P);
xlabel('距离(m)');
ylabel('幅值');
title(['声波传播 t=',num2str(t2)]);
end
end