数字信号处理

         实验2   傅里叶变换

1. 连续时间傅里叶变换

dt = 0.00005;                                                 

t = -0.005:dt:0.005;

A=100;a=50*sqrt(2)*pi;b=a;

     xa=exp(-a*t).*sin(b*t);                                     

%contunites_time fourier transform

Wmax = 2*pi* 2000;                                       

K = 500;                                                     

k = 0: 1: K;

W = k*Wmax/K;

Xa = xa * exp(-1i * t'*W)*dt;

Xa = real(Xa);

W = [-fliplr(W),W(2:501)];           %Flip matrix left to right  倒置 左右逐个交换

Xa = [fliplr(Xa),Xa(2:501)];            %xa over -Xa to Xa 合并矩阵

subplot(2,1,1);

plot(t*1000,xa);grid

title('analog signal');

xlabel('t  (ms)');

ylabel('xa(t)');

subplot(2,1,2);

plot(W/(2*pi*1000),Xa*1000);grid

title('continues_time fourier transform');

xlabel('f   (Khz)');

ylabel('Xa(jw) * 1000');

 

2. 离散时间傅里叶变换

dt = 0.00005;                             %时间步进量

t = -0.005:dt:0.005;  

A=100;a=50*sqrt(2)*pi;b=a;

xa=exp(-a*t).*sin(b*t);

ts = 0.0002;

n = -25:1:25;

x=exp(-a*n*ts).*sin(b*n*ts);

K = 500;

k = 0:1:K;

w = 2*pi*k/K;

X = x*exp(-1i *n'*w);

X= real(X);

w = [-fliplr(w),w(2:K+1)];

X = [fliplr(X),X(2:K+1)];                    %要对应-w所求的值

subplot(2,1,1);    

plot(t*1000,xa);

title('discrete signal');

xlabel('t in msec');

ylabel('x(n)');

hold on

stem(n*ts*1000,x);  

gtext('ts = 0.2msec'); hold off

subplot(2,1,2);

plot(w/pi, X);

title('discrete_time fourier transform');

xlabel('frequence in pi unit');

ylabel('X(w)');

 

3.快速傅里叶变换

法

Y = fft(X)

Y = fft(X,n)

Y = fft(X,n,dim)

Y = fft(X) 用快速傅里叶变换 (FFT) 算法计算 X 的离散傅里叶变换 (DFT)。

如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换。

如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。

如果 X 是一个多维数组，则 fft(X) 将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。

Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。

如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。

如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。

如果 X 是矩阵，则每列的处理与在向量情况下相同。

如果 X 为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。

Y = fft(X,n,dim) 返回沿维度 dim 的傅里叶变换。例如，如果 X 是矩阵，则 fft(X,n,2) 返回每行的 n 点傅里叶变换。

（1）计算信号的快速傅里叶变换

Fs=1000; %采样率

N = 1; % 周期数

t=N; % 信号时长 s

n=0:1/Fs:t-1/Fs; % 采样时间点，刚好采N 个周期

len = length(n); % 信号点数

y = sin(2*pi*20*n)+1.5*sin(2*pi*33*n)+0.6*sin(2*pi*50*n); % 采集到的离散信号

plot(n,y)

title('信号y')

 

Y = fft(y);

Y = abs(Y);

plot(Y(1:len/2));

title('直接调用fft')

 

clear all

f0 = 100;

fs = 500;

Ts = 1/fs;

n=1:1:100;

N = length(n);

y = sin(2*pi*f0*n*Ts);

plot(n,y)

y_fft=fft(y);

P2_y_fft =abs(y_fft/N);

P1_y_fft = P2_y_fft(1:N/2+1);

P1_y_fft(2:end-1) = 2*P1_y_fft(2:end-1)

f = fs*(0:N/2)/N;

figure

plot(f,P1_y_fft)

xlabel('f (Hz)')

ylabel('|P1(f)|')

 

（2） 使用傅里叶变换求噪声中隐藏的信号的频率分量

指定信号的参数，采样频率为 1 kHz，信号持续时间为 1.5 秒。

Fs = 1000;            % Sampling frequency                    

T = 1/Fs;             % Sampling period       

L = 1500;             % Length of signal

t = (0:L-1)*T;        % Time vector

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

X = S + 2*randn(size(t));

plot(1000*t(1:50),X(1:50))

title('Signal Corrupted with Zero-Mean Random Noise')

xlabel('t (milliseconds)')

ylabel('X(t)')

 

Y = fft(X);

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

plot(f,P1) 

title('Single-Sided Amplitude Spectrum of X(t)')

xlabel('f (Hz)')

ylabel('|P1(f)|')

 

Y = fft(S);

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1) 

title('Single-Sided Amplitude Spectrum of S(t)')

xlabel('f (Hz)')

ylabel('|P1(f)|')

 

1. 比较时域和频域中的余弦波

指定信号的参数，采样频率为 1kHz，信号持续时间为 1 秒。

Fs = 1000;                    % Sampling frequency

T = 1/Fs;                     % Sampling period

L = 1000;                     % Length of signal

t = (0:L-1)*T;                % Time vector

创建一个矩阵，其中每一行代表一个频率经过缩放的余弦波。结果 X 为 3×1000 矩阵。第一行的波频为 50，第二行的波频为 150，第三行的波频为 300。

x1 = cos(2*pi*50*t);          % First row wave

x2 = cos(2*pi*150*t);         % Second row wave

x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

在单个图窗中按顺序绘制 X 的每行的前 100 个项，并比较其频率。

for i = 1:3

    subplot(3,1,i)

    plot(t(1:100),X(i,1:100))

title(['Row ',num2str(i),' in the Time Domain'])end

 

出于算法性能的考虑，fft 允许您用尾随零填充输入。在这种情况下，用零填充 X 的每一行，以使每行的长度为比当前长度大的下一个最小的 2 的次幂值。使用 nextpow2 函数定义新长度。

n = 2^nextpow2(L);

指定 dim 参数沿 X 的行（即对每个信号）使用 fft。

dim = 2;

计算信号的傅里叶变换。

Y = fft(X,n,dim);

计算每个信号的双侧频谱和单侧频谱。

P2 = abs(Y/L);

P1 = P2(:,1:n/2+1);

P1(:,2:end-1) = 2*P1(:,2:end-1);

在频域内，为单个图窗中的每一行绘制单侧幅值频谱。

for i=1:3

    subplot(3,1,i)

    plot(0:(Fs/n):(Fs/2-Fs/n),P1(i,1:n/2))

title(['Row ',num2str(i),' in the Frequency Domain'])

end