Table of contents
EditEditEdit%7.2.4 Verify the scaling properties of the Fourier transform
%7.3.2 Effects of amplitude distortion on hearing and vision
%7.3.3 Effects of phase distortion on hearing and vision
[Purpose]
1. Learn to use MATLAB to plot the logarithmic amplitude-frequency characteristics and phase-frequency characteristics of the frequency response function.
2. Learn to use MATLAB to complete signal sampling and spectral analysis of sampled signals.
3. Learn to use MATLAB to reconstruct the sampled signal.
4. Understand the application of other Fourier analysis using MATLAB.
[experiment apparatus]
- computer
- MATLAB software
[Experiment content]
1. The frequency response function of a system
, try to draw its logarithmic amplitude-frequency characteristics and phase-frequency characteristics.
omega=logspace(-3,1,500); %logspace function to generate log-spaced vectors
h=1./(1+1i*2*omega);
figure;
subplot(2,1,1);grid on;
semilogx(omega,20*log10(abs(h)));%exist x axis logarithm
title('Logarithmic Amplitude Frequency Characteristics');
subplot(2,1,2);grid on;
semilogx(omega,angle(h));
title('Phase frequency characteristics');
2. Try to draw the frequency response function
The logarithmic amplitude-frequency characteristic of .
w=0:0.02:100;
magHw=abs(3./(3*1i*w+2-w.^2));
semilogx(w,magHw);
title('Logarithmic Amplitude Frequency Characteristics')
xlabel('Frequency in rad/sec-log scale');
ylabel('Magnitude of Vout/Vin');
grid
3. The known signal is
, use MATLAB programming to realize the sampling signal fs(t) and its spectrum obtained by the impulse pulse of the signal. Let the parameters E=5, τ=0.5, use the sampling interval
s=1;
dt=0.1;
Ts=0.2;
t1=-4:dt:4;
ft=(5*(1+cos(2*pi*t1))/2).*(heaviside(t1+0.5)-heaviside(t1-0.5));%Expression of step function
subplot(221);plot(t1,ft),
grid on;
axis([-2 2 -0.1 5.1]);
xlabel('Time(sec)'),ylabel('f(t)');
title('Signal ft');
N=500;
k=-N:N;
W=pi*k/(N*dt);
Fw=dt*ft*exp(-1i*t1'*W);
subplot(222);plot(W,abs(Fw)),
grid on;axis([-30 30 -0.1 1.1*pi]);
xlabel('\omerga'),ylabel('F(w)');
title('Signal ft spectrum');
Ts=0.2;
t2=-4:Ts:4;
fst=(5*(1+cos(2*pi*t2))/2).*(heaviside(t2+0.5)-heaviside(t2-0.5));
subplot(223);plot(t1,ft,':'),hold on;
stem(t2,fst),grid on;
axis([-2 2 -0.1 5.1]);
xlabel('Time(sec)'),ylabel('fs(t)');
title('sampled signal'),hold off;
Fsw=Ts*fst*exp(-1i*t2'*W);
subplot(224);plot(W,abs(Fsw)),grid on;
axis([-30 30 -0.2 1.1*pi]);
xlabel('\omega'),ylabel('Fs(w)');
title('Spectrum of the sampled signal')
4. For the sampled signal obtained in question 3, use a low-pass filter with a cutoff frequency of 4pi to filter it to reconstruct the signal f(t), and calculate the absolute error between the reconstructed signal and the original raised cosine pulse signal.
E=5;
wc=4*pi;
Ts=0.2;
n=-100:100;
nTs=n*Ts;
fs=(E*(1+cos(2*pi*nTs))/2).*(heaviside(nTs+0.5)-heaviside(nTs-0.5));
t=-4:0.1:4;
ft=fs*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
t1=-4:0.1:4;
f1=(5*(1+cos(2*pi*t1))/2).*(heaviside(t1+0.5)-heaviside(t1-0.5));
subplot(311);
plot(t1,ft,':'),hold on;
stem(nTs,fs),grid on;
axis([-2 2 -0.1 5]);
xlabel('nTs'),ylabel('f(nTs)');
title('sampling interval Ts sampled signal f(nTs)');
hold off;
subplot(312);
plot(t,ft),grid on;
axis([-4 4 -0.1 5]);
xlabel('t'),ylabel('f(t)');
title('Depend on f(nTs)Signal reconstruction to get the original signal');
error=abs(ft-f1);subplot(313);
plot(t,error),grid on;
xlabel('t'),ylabel('error(t)');
title('The absolute error between the reconstructed signal and the original signal')
5. Knowing the modulation signal f(t)=cos5πt and the carrier signal fc(t)=cos60πt, program to draw the waveform diagram and spectrogram during the modulation and demodulation process.
clear clc
fs=40;
Fs=400;
Fc=30;
N=400;
n=0:N;
t=n/Fs;
xt=cos(pi*5*t);
xct=cos(2*pi*Fc*t);
yt=xt.*xct;
Xw=fftshift(abs(fft(xt,512)));
Yw=fftshift(abs(fft(yt,512)));
ww=-256:255;
ww=ww*Fs/512;
figure,subplot(2,1,1),plot(t,xt),title('Modulated Signal Waveform');
subplot(2,1,2),plot(ww,Xw),title('Modulated Signal Spectrum');
figure,subplot(2,1,1),plot(t,yt),title('Modulated Signal Waveform');
subplot(2,1,2),plot(ww,Yw),title('Modulated Signal Spectrum');
ylt=yt.*xct;
figure,subplot(2,1,1),plot(t,ylt),title('Intermediate signal waveform during demodulation');
Ylw=fftshift(abs(fft(ylt,512)));
subplot(2,1,2),plot(ww,Ylw),title('Intermediate signal spectrum during demodulation');
h=fir1(20,40/200,hamming(21));
figure,freqz(h,1),title('Filter Frequency Characteristics');
y2t=filter(h,1,ylt);
Y2w=fftshift(abs(fft(y2t,512)));
figure,subplot(2,1,1),plot(t,y2t),title('The waveform of the demodulated signal');
subplot(2,1,2),plot(ww,Y2w),title('The spectrum of the demodulated signal')
%Run the program, observe the results, and carefully complete the following questions.
% ①How to understand that any periodic signal can be decomposed into the form of the weighted sum of sinusoidal signals;
% Answer: Any continuous measurement sequence or signal can be expressed as an infinite superposition of sine wave signals of different frequencies. The Fourier series is the weight of the weight. The Fourier series decomposition method can express any periodic signal as the sum of the sine and cosine signals, and the cosine can be expressed as a sine with a phase shift of 90°. It can also be clearly seen through experiments that when the maximum harmonic order is 400, the signal tends to be a periodic square wave signal.
% ②How to understand that the Fourier transform of aperiodic signal can be obtained by normalizing the Fourier coefficient of the periodic signal by its fundamental frequency;
% Answer: The Fourier transform of aperiodic signals can be regarded as a periodic signal with an infinite period. At this time, the spectral interval 2*pi/w tends to zero, and the spectral lines become continuous. In order to compare the ratios between infinitesimals, the concept of normalization is introduced, and the fundamental frequency of the Fourier coefficients is normalized in the image. After normalization, the frequency domain energy image can be better quantized (re- Grayscale images quantized to the 0-255 range)
% ③ Reasonably explain that the high-frequency information of the signal corresponds to its details, and the low-frequency information of the signal corresponds to its overview;
% Answer: A two-dimensional image can be decomposed into different frequency components. Among them, low-frequency components describe a wide range of information, while high-frequency components describe specific details. In a grayscale image, areas with small changes in brightness are mainly low-frequency components, while areas with sharp changes in brightness (such as the edges of objects) are mainly high-frequency components. The low-frequency signal corresponds to the slow-changing part of the information, and the high-frequency signal corresponds to the slow-changing part of the information. The general outline changes slowly, and the details change faster. Therefore, the high-frequency information of the signal corresponds to its details. The low-frequency information corresponds to its profile.
% ④ Explain the reason why the time domain shift signal does not affect the amplitude spectrum of its frequency domain;
% Answer: If the signal is shifted left/right by t0 in the time domain, the frequency domain signal will be multiplied by exp(±jwt0), and its amplitude is 1, that is, its amplitude spectrum will not be affected.
% ⑤ Describe the impact of audio signals on hearing after scale transformation in the time domain, and explain the reasons;
% Answer: After the audio signal is scale-expanded in the time domain, the sound becomes gentle, and after the scale is compressed, the sound becomes short. This affects the loudness of the audio, in the form of lower volume after expansion and higher volume after compression
% ⑥According to the spectrogram of the audio signal in the time domain differential signal, explain the reason why the sound heard is different from the original sound (before differentiation);
% Answer: After differentiating the signal, the frequency domain signal is equivalent to multiplying by jw, the amplitude of the high frequency signal increases, and the amplitude of the low frequency signal decreases. Therefore, the differentiated signal has more high-frequency components and less low-frequency components, which has a certain sharpening effect.
% ⑦After the sound signal is shifted by 1Hz in the frequency domain, according to the time domain diagram corresponding to the inverse Fourier transform modulo signal, explain the change of the heard sound, and explain the reason in theory;
% Answer: The size of the sound will fluctuate to a certain extent, because the frequency domain shift is equivalent to multiplying the time domain signal by exp(jwt), which is multiplied by the trigonometric function in the real number domain, and its size will fluctuate periodically.
% ⑧ Discuss the effect of phase distortion on human hearing and vision;
% Answer: Visually, the edge of the two-dimensional image seen is less sharp, the image is more blurred, the difference between different gray values is also reduced, and the auditory performance is that the difference between different frequencies of sound becomes smaller.
% ⑨ After the sound signal is reconstructed into a two-dimensional image signal, discuss the difference between musical sound and noise on visual effects;
% Answer: There are many high-frequency signals in the music signal, and the corresponding image details are more, and the low-frequency signals of noise are more, and the image is smoother.
% ⑩ After the image signal is expanded into a one-dimensional signal, the differences in auditory effects of one-dimensional audio signals corresponding to regular texture images, arbitrary images and noise images are discussed.
% Answer: There are many regular parts in the audio signal converted from regular texture images, which may form some fixed repeating parts such as harmonics and rotations, which makes people impressive, while the randomness of any image increases, and the sense of regularity increases. will be reduced, and the noise image has a poor auditory effect due to more cluttered parts.
% 7.1.2 Using the first 4, 40 and 400 items respectively, draw an approximate diagram of the periodic rectangular pulse signal
clear;close all; format compact
t=-1:0.001:1;
E=1;T=1;w=2*pi/T;
ft=square(2*pi*t,50); %Write a function expression
figure(1);
h=figure(1);
set(h,'name','7.1.2 Using the first 4, 40, and 400 items respectively, draw an approximate diagram of the periodic rectangular pulse signal');
set(h,'position', [0.2 0.3 0.6 0.6],'Units','normalized' );%Plot in relative size
set(h,'position', [0.2 0.3 0.6 0.6],'Units','normalized' );
subplot(2,2,1),plot(t,ft),grid on
axis([-1,1,-1.5,1.5])%set axis
title('Periodic square wave signal')
n_max=[4 40 400];
N=length(n_max);
for k=1:N %pass for Drawing in a loop
n=1:2:n_max(k);
b=4*E./(pi*n); %set amplitude
x=b*sin(2*pi*n'*t); %set function expression
subplot(2,2,k+1),plot(t,x),grid on
axis([-1,1,-1.5,1.5])
title(['Maximum harmonic order=',num2str(n_max(k))])
end %end loop flag
pause;
%7.1.4 Time Progressive Synthesis Demonstration of Calculation and Reconstruction Algorithm of Fourier Coefficients of Periodic Kernel Function of Arbitrary Periodic Signals
Um=0.1;T=1;w=2*pi/T; %Initialize max, period and frequency separately
num_noise=100; %define noise
%N=input('Take the harmonic order N= ');
N=60; %Before setting N item
Fs=1000;
t=0:1/Fs:T;
num_points=numel(t); %Set the number of time points
dt=T/(num_points-1); %Set time point interval
un=Um*(sin(1*2*pi*t)); %Decomposed waveform
n=randn(num_noise-2,1); %Randomly generated matrix
n=[0;n;0];
n = interp1([0:num_noise-1],n,linspace(0,num_noise-2,num_points),'linear'); %Linear interpolation is used here to insert the previous noise signal
un=un+0.8*n;
figure(2);
h=figure(2);
set(h,'position', [0.2 0.05 0.6 0.8],'Units','normalized' ); %Drawing in relative size
set(h,'position', [0.2 0.05 0.6 0.8],'Units','normalized' );
subplot(3,1,1);plot(t,un,'LineWidth',3);
axis([-0.05*T T+0.05*T -1.8 1.8]);
grid on;
hold on;
a0=trapz(un)*dt/T*1; %Initialization parameters a0,b0,Phi0
b0=0;
A0=abs(a0);
Phi0=0;
for k=1:N %Solve in a loop a(k) ,b(k), A(k), Phi(k)
a(k)=trapz(un.*cos(k*w*t))*dt/T*2;
b(k)=trapz(un.*sin(k*w*t))*dt/T*2;
A(k)=sqrt(a(k)^2+b(k)^2);
Phi(k)=atan2(-b(k),(a(k)+eps));
end
subplot(3,1,2);stem(0,A0); axis([0 N -inf,inf]);title('One-sided amplitude spectrum')
subplot(3,1,3);stem(0,Phi0); axis([0 N -inf,inf]);title('one-sided phase spectrum')
Reu=A0;
subplot(3,1,1);plot(t,Reu+0*t,'g');
for k=1:N
clf;
subplot(3,1,1);plot(t,un,'LineWidth',3);grid on;hold on;
axis([-0.01*T T+0.01*T -1.8 1.8]); %Initialize the axes
title(strcat('original and before',num2str(k),'Sum of Subharmonics'));
subplot(3,1,2);stem(0:k,[A0,A(1:k)]);title('One-sided amplitude spectrum');axis([0 N -inf,inf]);grid on;
subplot(3,1,3);stem(0:k,[Phi0,Phi(1:k)]);title('one-sided phase spectrum');axis([0 N -3.2 3.2]);grid on;
ad=A(k)*cos(k*w*t+Phi(k));
subplot(3,1,1);plot(t,ad,'k');
subplot(3,1,1);plot(t,Reu,'r','LineWidth',2);
if k==18
pause;
end
pause(0.1);
Reu=Reu+ad;
end
Psll=trapz(un.^2)*dt/T; % original signal energy
Ps12=A0^2+sum(A(1:end).^2/2); % forward N term harmonic energy
pause;
%7.2.2 Compare the relationship between Fourier transform of aperiodic signal and Fourier coefficient of periodic signal
figure(3);
h=figure(3);
set(h,'name','7.1.4 Calculation of Fourier Coefficients of Periodic Kernel Function of Arbitrary Periodic Signals and Demonstration of Time Progressive Synthesis of Reconstruction Algorithms');
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' ); %Drawing in relative size
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
subplot(2,1,1);stem(0:N,[A0*2,A(1:end)]); title('One-sided amplitude spectrum');grid on;
subplot(2,1,2);stem(0:N,[Phi0,Phi(1:end)]);title('one-sided phase spectrum');grid on;
uf=fft(un,numel(un))/(Fs*T); %Because each cycle has Fs points
figure(4);
h=figure(4);
set(h,'position', [0.2 0.5 0.6 0.4],'Units','normalized' ); %Drawing in relative size
set(h,'position', [0.2 0.5 0.6 0.4],'Units','normalized' );
set(h,'name','7.2.2 Comparing the relationship between Fourier transform of aperiodic signal and Fourier coefficient of periodic signal');
subplot(2,1,1);stem(0:N,abs(uf(1:N+1)));title('Aperiodic signal amplitude spectrum (positive frequency part)');grid on;
subplot(2,1,2);stem(0:N,angle(uf(1:N+1)));title('Aperiodic signal phase spectrum (positive frequency part)');grid on;
pause;
%7.2.3 Use the MATLAB symbolic operation function fourier to solve the Fourier transform of the following signals, and use ifourier to verify
syms t;
ft1=exp(-t)*heaviside(t);%set function
ft2=1/(sqrt(2*pi))*exp(-((t^2)/2));%set function
Fw1=fourier(ft1);%Fourier transform of a function
ft1i=ifourier(Fw1);%Inverse Fourier Transform of a Function
Fw2=fourier(ft2);
ft2i=ifourier(Fw2);
figure(7);
h=figure(7);
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
set(h,'position', [0.2 0.05 0.6 0.4],'Units','normalized' );
set(h,'name','7.2.3 use MATLAB Symbolic operation functions fourier Solve the Fourier transform of the following signal using ifourier verify');
subplot(2,1,1);ezplot(abs(Fw1));title('Amplitude spectrum of original signal 1');axis([-6 6 0 1.1]); grid on;
subplot(2,1,2);h1=ezplot(abs(Fw2),[-6 6]);title('Amplitude spectrum of original signal 2');axis([-6 6 0 1.1]); grid on;
set(h1,'Color','r','LineWidth',2) %Set the description line to red
figure(5);
h=figure(5);
set(h,'position', [0 0.5 0.5 0.4],'Units','normalized' );
set(h,'position', [0 0.5 0.5 0.4],'Units','normalized' );
set(h,'name','7.2.3 use MATLAB Symbolic operation functions fourier Solve the Fourier transform of the following signal using ifourier verify');
subplot(2,1,1);ezplot(ft1); title('original time domain signal');axis([-0.1 3 0 1.1]); grid on;
subplot(2,1,2);h1=ezplot(ft1i); title('Inversely transformed time domain signal');axis([-0.1 3 0 1.1]); grid on;%Image Verification via Inverse Fourier Transform
set(h1,'Color','c','LineWidth',2);%Set the description line to blue
figure(6);
h=figure(6);
set(h,'position', [0.5 0.5 0.5 0.4],'Units','normalized' );
set(h,'position', [0.5 0.5 0.5 0.4],'Units','normalized' );
set(h,'name','7.2.3 use MATLAB Symbolic operation functions fourier Solve the Fourier transform of the following signal using ifourier verify');
subplot(2,1,1);h1=ezplot(ft2);title('original time domain signal');axis([-3 3 0 0.5]); grid on;
set(h1,'Color','r');
subplot(2,1,2);h1=ezplot(ft2i);title('Inversely transformed time domain signal');axis([-3 3 0 0.5]); grid on;
set(h1,'Color','c','LineWidth',2);%Image Verification via Inverse Fourier Transform
pause;
%7.2.4 Verify the scaling properties of the Fourier transform
dt=0.1;
t1=-4:dt:4;
ft1=((1+cos(t1))/2).*(heaviside(t1+pi)-heaviside(t1-pi));%define functional ft1
ft2=((1+cos(2*t1))/2).*(heaviside(t1+0.5*pi)-heaviside(t1-0.5*pi));%define functional ft2
N=500;
k=-N:N;
W=pi*k/(N*dt);
Fw1=dt*ft1*exp(-1i*t1'*W);%respectively ft1 and ft2 take a Fourier transform
Fw2=dt*ft2*exp(-1i*t1'*W);
figure(8);
h=figure(8);
set(h,'position', [0 0.1 0.5 0.8],'Units','normalized' );
set(h,'position', [0 0.1 0.5 0.8],'Units','normalized' );
set(h,'name','7.2.4 Verify the scaling properties of the Fourier transform');
subplot(2,1,1);plot(t1,ft1,'LineWidth',2); title('Signal time domain waveform'); axis([-4 4 -0.1 1.1]); grid on;%Plot signal time domain waveform
subplot(2,1,2);plot(t1,ft2,'r','LineWidth',2); title('Scaled time-domain waveform');axis([-4 4 -0.1 1.1]); grid on;%Plot the scaled signal time-domain waveform
figure(9);
h=figure(9);
set(h,'position', [0.5 0.1 0.5 0.8],'Units','normalized' );
set(h,'position', [0.5 0.1 0.5 0.8],'Units','normalized' );
set(h,'name','7.2.4 Verify the scaling properties of the Fourier transform');
subplot(2,1,1);plot(W,real(Fw1),'LineWidth',2);title('composite spectrum of the signal');axis([-8 8 -0.6 3.5]); grid on;%Plot the composite spectrum of a signal
subplot(2,1,2);plot(W,real(Fw2),'r','LineWidth',2);title('Composite spectrum of the scaled signal');axis([-8 8 -0.6 3.5]); grid on;%Plot the scaled signal composite spectrum
pause;
%7.3.2 Effects of amplitude distortion on hearing and vision
T=3;
Fs=40000;
t=0:1/Fs:T;
s1=sin(600*2*pi*t);%set separately s1，s2，s3
s2=sin(850*2*pi*t);
s3=3*sin(850*2*pi*t);
s1s2=s1+s2;%Define the form of the sum of functions
s1s3=s1+s3;
figure(11);
h=figure(11);
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'name','7.3.2 Effects of Amplitude Distortion on Hearing and Vision');%draw s1,s2,s3,s1s2,s1s3 time domain diagram of
subplot(511);plot(t(1:600),s1(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1');
subplot(512);plot(t(1:600),s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s2');
subplot(513);plot(t(1:600),s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s3');
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
[R,C]=meshgrid(0:1/Fs:600/160000:1/Fs:600/160000);
s1i=sin(600*2*pi*R);%set separately s1i，s2i，s3i
s2i=sin(850*2*pi*R);
s3i=3*sin(850*2*pi*R);
s1s2i=sin(600*2*pi*R)+sin(850*2*pi*R);
s1s3i=sin(600*2*pi*R)+3*sin(850*2*pi*R);
figure(13);
h=figure(13);
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'name','7.3.2 Effects of Amplitude Distortion on Hearing and Vision');
subplot(5,1,1),imshow(s1i(1:30,:)/3,[0,1]);title('s1i')
subplot(5,1,2),imshow(s2i(1:30,:)/3,[0,1]);title('s2i')
subplot(5,1,3),imshow(s3i(1:30,:),[]);title('s3i')
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i');
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i');
pause;
figure(11);
subplot(514);plot(t(1:600),s1s2(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(13); %Display the signal in the form of sound
subplot(5,1,4),imshow(mat2gray(s1s2i(1:30,:))*64,spring);title('s1s2i');
sound(s1s2,Fs);
pause(4);
figure(11);
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(13);
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i');
figure(11);
subplot(515);plot(t(1:600),s1s3(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(13);
subplot(5,1,5),imshow(mat2gray(s1s3i(1:30,:))*64,spring);title('s1s3i');
sound(s1s3,Fs);
pause(4);
figure(11);
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(13);
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i');
pause;
%7.3.3 Effects of phase distortion on hearing and vision
T=3;
Fs=40000;
t=0:1/Fs:T;
s1=sin(600*2*pi*t);%set separately s1，s2，s3 and the sum of functions
s2=sin(850*2*pi*t);
s3=sin(850*2*pi*t-0.7);
s1s2=s1+s2;
s1s3=s1+s3;
figure(14);
h=figure(14);
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.6 0.85],'Units','normalized' );
set(h,'name','7.3.3 Effects of phase distortion on hearing and vision');
subplot(511);plot(t(1:600),s1(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1');
subplot(512);plot(t(1:600),s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s2');
subplot(513);plot(t(1:600),s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s3');
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
[R,~]=meshgrid(0:1/Fs:600/160000:1/Fs:600/160000);
s1i=sin(600*2*pi*R);%set separately s1i，s2i，s3i and the sum of functions
s2i=sin(850*2*pi*R);
s3i=sin(850*2*pi*R-0.7);
s1s2i=sin(600*2*pi*R)+sin(850*2*pi*R);
s1s3i=sin(600*2*pi*R)+3*sin(850*2*pi*R-0.7);
figure(16);
h=figure(16);
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'position', [0.62 0.05 0.37 0.85],'Units','normalized' );
set(h,'name','7.3.3 Effects of phase distortion on hearing and vision');
subplot(5,1,1),imshow(s1i(1:30,:),[0,1]);title('s1i')%draw images separately
subplot(5,1,2),imshow(s2i(1:30,:),[0,1]);title('s2i')
subplot(5,1,3),imshow(s3i(1:30,:),[]);title('s3i')
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i');
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i');
pause;
figure(14);
subplot(514);plot(t(1:600),s1s2(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(16);
subplot(5,1,4),imshow(mat2gray(s1s2i(1:30,:))*64,spring);title('s1s2i');
sound(s1s2,Fs);
pause(4);
figure(14);
subplot(514);plot(t(1:600),s1s2(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s2');
figure(16);
subplot(5,1,4),imshow(s1s2i(1:30,:),[]);title('s1s2i');
figure(14);
subplot(515);plot(t(1:600),s1s3(1:600),'r');axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(16);
subplot(5,1,5),imshow(mat2gray(s1s3i(1:30,:))*64,spring);title('s1s3i');
sound(s1s3,Fs);
pause(4);
figure(14);
subplot(515);plot(t(1:600),s1s3(1:600));axis([0,600/Fs,-3.7,3.7]);title('s1s3');
figure(16);
subplot(5,1,5),imshow(s1s3i(1:30,:),[]);title('s1s3i');
pause;
%7.4.2 Time-shift the audio signal and compare it with the original signal in the time domain and frequency domain
Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%read in audio
testsou1=testsou1(:,1);%Split the audio
testsou1=resample(testsou1,Fs,Fs1);%change its dimensions
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A2=[zeros(Fs,1);A1];%time shift one second
A=[A1;zeros(Fs,1)];
Af=fft(A,Fs*5); %Using the Fast Fourier Convolution Function
A2f=fft(A2,Fs*5);
figure(18);
h=figure(18);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.2 Time-shift the audio signal and compare it with the original signal in the time and frequency domains');
subplot(2,1,1), plot([1:round(Fs*5/2)],abs(Af(1:round(Fs*5/2)))); title('Af')% draw its image
subplot(2,1,2), plot([1:round(Fs*5/2)],abs(A2f(1:round(Fs*5/2)))); title('A2f')
figure(17);
h=figure(17);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.2 Time-shift the audio signal and compare it with the original signal in the time and frequency domains');
subplot(2,1,1), plot(A), axis([0,numel(A), -1, 1]);grid on;title('A')% draw its image
subplot(2,1,2), plot(A2), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause;
sound([A,A2],Fs);
subplot(2,1,1), plot(A,'r'), axis([0,numel(A), -1, 1]);grid on;title('A')% draw its image
subplot(2,1,2), plot(A2,'r'), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause(6.6);
subplot(2,1,1), plot(A), axis([0,numel(A), -1, 1]);grid on;title('A')% draw its image
subplot(2,1,2), plot(A2), axis([0,numel(A2), -1, 1]);grid on;title('A2')
pause;
%7.4.3 Scale the audio signal and compare it with the original signal in the time domain and frequency domain
Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%read audio
testsou1=testsou1(:,1);%Split the audio
testsou1=resample(testsou1,Fs,Fs1);%Change the sampling rate to implement decimation and interpolation operations
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A3=resample(A1,numel(A1)*0.6,numel(A1));%scaling a=1/0.6
A1f=fft(A1,Fs*5); %Use the Fast Fourier Transform function
A3f=fft(A3,Fs*5);
figure(20);
h=figure(20);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.3 Scale the audio signal and compare it with the original signal in the time and frequency domains');
subplot(2,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')%draw the resulting image
subplot(2,1,2), plot([1:round(Fs*5/2)],abs(A3f(1:round(Fs*5/2)))); title('A3f')
figure(19);
h=figure(19);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.3 Scale the audio signal and compare it with the original signal in the time and frequency domains');
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1') %draw the resulting image
subplot(2,1,2), plot(A3), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause;
sound(A1,Fs);
subplot(2,1,1), plot(A1,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A1')
pause(5)
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1') %draw the resulting image
pause(1.5)
sound(A3,Fs);
subplot(2,1,2), plot(A3,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause(4)
subplot(2,1,2), plot(A3), axis([0,numel(A1), -1, 1]);grid on;title('A3')
pause;
%7.4.4 Perform primary and secondary differential operations on the audio signal respectively, and compare the amplitude spectrum and sound changes of the two signals with the original signal.
Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%read in audio
testsou1=testsou1(:,1);%Split the audio
testsou1=resample(testsou1,Fs,Fs1);%to sample
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
A4=[0;A1(1:end-1)]-A1; %find first derivative
A5=[0;A4(1:end-1)]-A4; %Find the second derivative
A1f=fft(A1,Fs*5);
A4f=fft(A4,Fs*5);
A5f=fft(A5,Fs*5);
figure(21);
h=figure(21);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','Perform primary and secondary differential operations on audio signals respectively');
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')%Perform first-order and second-order differentiation operations on audio signals and draw images
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2)))); title('A4f')
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2)))); title('A5f')
pause;
sound(A1,Fs);
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2))),'r'); title('A1f')
pause(6)
subplot(3,1,1), plot([1:round(Fs*5/2)],abs(A1f(1:round(Fs*5/2)))); title('A1f')
pause(0.5)
sound(A4,Fs);
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2))),'r'); title('A4f')
pause(6)
subplot(3,1,2), plot([1:round(Fs*5/2)],abs(A4f(1:round(Fs*5/2)))); title('A4f')
pause(0.5)
sound(A5,Fs);
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2))),'r'); title('A5f')
pause(6)
subplot(3,1,3), plot([1:round(Fs*5/2)],abs(A5f(1:round(Fs*5/2)))); title('A5f')
pause;
%7.4.5 Perform a frequency-shift 1Hz operation on the Fourier transform of the audio signal, and compare the time-domain waveform and sound changes
Fs=8000;
[testsou1,Fs1]=audioread('aud1.mp3');%read in audio
testsou1=testsou1(:,1);%Split the audio
testsou1=resample(testsou1,Fs,Fs1);%sampling
tem=testsou1(120000:120000+5*8000-1);
A1=tem;
Af=fft(A1);%Use the Fast Fourier Transform function
Afsh=[zeros(round(numel(Af)/Fs),1);Af(1:round(numel(Af)/2));...
Af(round(numel(Af)/2)+1:end);zeros(round(numel(Af)/Fs),1)];%Frequency Domain 1 Hz
Ash=ifft(Afsh,40000);%Frequency-shifted time-domain signal, inversely transformed
figure(22);
h=figure(22);
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.01 0.05 0.48 0.85],'Units','normalized' );
set(h,'name','7.4.5 Frequency-shifting the Fourier transform of an audio signal by 1 Hz operate');
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1')
subplot(2,1,2), plot(real(Ash)), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
figure(23);
h=figure(23);
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
set(h,'position', [0.51 0.05 0.48 0.85],'Units','normalized' );
plot([200:450],abs(Af(200:450)),'k');
hold on
plot([200:450],abs(Afsh(200:450)),'m','LineWidth',2);
grid on;
title('Af and Afsh')
pause;
sound(A1,Fs);
figure(22);
subplot(2,1,1), plot(A1,'r'), axis([0,numel(A1), -1, 1]);grid on;title('A1')%draw its image
pause(5)
subplot(2,1,1), plot(A1), axis([0,numel(A1), -1, 1]);grid on;title('A1')
pause(1)
sound(real(Ash),Fs);
subplot(2,1,2), plot(real(Ash),'r'), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
pause(5)
subplot(2,1,2), plot(real(Ash)), axis([0,numel(A1), -1, 1]);grid on;title('Ash')
pause; 
%7.5 sees sound
[testsou,Fs1]=audioread('aud1.mp3');%Read in audio file
testsou1=testsou(:,1);%Split the audio
Ig1=reshape(testsou1(1:100000),[250,400]); %Reorganize the resulting matrix into a new matrix
Ig1f=fft(testsou1,10000);%Fast Fourier Transform it
noise1=randn(100000,1);%generate noise
ng1=reshape(noise1,[250,400]);%Reorganize the noise into a new matrix
ng1f=fft(noise1,10000);%Fast Fourier Transform it
figure(24);
h=figure(24);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );%Plot audio grayscale values
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','7.5 see the sound');
subplot(2,2,1),imshow(Ig1,[]);title('Ig1')
subplot(2,2,2),imshow(ng1,[]);title('ng1')
subplot(2,2,3),plot(abs(Ig1f(1:5000)));title('Ig1f')
subplot(2,2,4),plot(abs(ng1f(1:5000)));title('ng1f')
pause;
sound([Ig1(:);zeros(40000,1);0.2*ng1(1:size(ng1,2)*size(ng1,1)*0.75)'],Fs1); %Play tones and noise in sequence
subplot(2,2,1),imshow(mat2gray(Ig1)*64,summer);title('Ig1');
pause(5);
subplot(2,2,1),imshow(Ig1,[]);title('Ig1');
pause(1);
sound(0.2*ng1(1:size(ng1,2)*size(ng1,1)*0.75),Fs1); %Play tones and noise in sequence
subplot(2,2,2),imshow(mat2gray(ng1)*64,summer);title('ng1');
pause(4);
subplot(2,2,2),imshow(ng1,[]);title('ng1');
pause;
%7.6 Hear the image
function g=gscale(f,varargin) %right gscale function to define
if length(varargin)==0
method='full8';
else method=varargin{1};
end
if strcmp(class(f),'double') & (max(f(:))>1 || min(f(:))<0)
f=mat2gray(f);
end
switch method
case 'full8'
g=im2uint8(mat2gray(double(f)));
case 'full16'
g=im2uint16(mat2gray(double(f)));
case 'minmax'
low = varargin{2};
high = varargin{3};
if low>1 || low<0 || high>1 || high<0
error('Parameters low and high must be in the range [0,1]')
end
if strcmp(class(f),'double')
low_in=min(f(:));
high_in=max(f(:));
elseif strcmp(class(f),'uint8')
low_in=double(min(f(:)))./255;
high_in=double(max(f(:)))./255;
elseif strcmp(class(f),'uint16')
low_in=double(min(f(:)))./65535;
high_in=double(max(f(:)))./65535;
end
g=imadjust(f,[low_in high_in],[low high]);
otherwise
error('Unknown method')
end
Fs=8000;
I1=imread('DSC_0333.JPG');%read in picture
%I2=imread('texture 7.JPG');
[R, C]=meshgrid(-2:0.01:2,-2:0.01:2);
u0=5*pi; v0=2*pi; A=1;
f1 = A*sin(u0*R + v0*C);%respectively f1，f2，f3 define
f2 = A*sin(u0*R) + A*sin(v0*C);
f3 = A*sin(u0*sqrt(R.^2+C.^2));
Ig2=gscale(f2); %Scale the input image using the function defined above
Ig1=rgb2gray(I1);%Will RGB Convert image to grayscale
Ig1=mat2gray(Ig1);%Normalization to Image Matrix
Ig1=imresize(Ig1,[400,800]);%Modify size
%Ig2=rgb2gray(I2);
Ig2=mat2gray(Ig2);
Ig2=imresize(Ig2,[400,800]);
Ig2fil=imfilter(Ig2,ones(15)/225); %filter operation
Ig3=randn(400,800);
figure(25),
h=figure(25);
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'position', [0.1 0.05 0.8 0.85],'Units','normalized' );
set(h,'name','7.6 hear the image');
subplot(3,1,1),imshow(Ig1,[]);title('Ig1')
subplot(3,1,2),imshow(Ig2,[]);title('Ig2')
subplot(3,1,3),imshow(Ig3,[]);title('Ig3')
pause;
sound([Ig1(:);zeros(40000,1);Ig2fil(:);zeros(40000,1);0.2*Ig3(:)],Fs*8);
subplot(3,1,1),imshow(mat2gray(Ig1)*64,summer);title('Ig1')%draw its image
pause(5.6);
subplot(3,1,1),imshow(Ig1,[]);title('Ig1')
subplot(3,1,2),imshow(mat2gray(Ig2)*64,summer);title('Ig2')
pause(5.4);
subplot(3,1,2),imshow(Ig2,[]);title('Ig2')
subplot(3,1,3),imshow(mat2gray(Ig3)*64,summer);title('Ig3')
pause(5.5);
subplot(3,1,3),imshow(Ig3,[]);title('Ig3')
[Experimental perception]
Through this experiment, I learned to plot the logarithmic amplitude-frequency and phase-frequency characteristics of the function, sample the signal, analyze the spectrum of the sampled signal, and reconstruct the sampled signal. And through experiments, I deepened my understanding of the course of signal processing, and learned the characteristics of signals through visual methods, including audio signals and image signals. The reconstruction of the signal allows me to see the difference between musical tones and noise, and through experiments, I have a more vivid understanding that any periodic signal can be decomposed into the form of sinusoidal signal weighting. And this experiment also gave me a deeper understanding of high-frequency and low-frequency signals from different dimensions (audio, image, etc.). In addition, I also learned how to use functions related to signal processing by searching documents. For example, imshow(I,[low high]) function is used to display grayscale image I, and specify the display range in the form of two element vector [low high], y = resample (x, P, q) uses a polyphase filter to resample the sequence in the vector x at p/q times the original sampling rate, the reshape function reshapes the original matrix into a new matrix, fft Fast Fourier convolution functions and so on, which also gave me a benefit from my engineering skills.