Используйте преобразования Фурье, чтобы найти частотные составляющие сигнала, зарытого в шум.

Задайте параметры сигнала с частотой дискретизации 1 кГц и длительностью сигнала 1,5 секунды.

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

Сформировать сигнал, содержащий синусоиду 50 Гц амплитуды 0,7 и синусоиду 120 Гц амплитуды 1.

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

Повреждите сигнал с нулевым средним белым шумом с отклонением 4.

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

Постройте график сигнала с шумом в временной интервал. Трудно идентифицировать частотные составляющие, посмотрев на сигнал X(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. Затем вычислите односторонний спектр P1 на основе P2 и четную длину сигнала L.

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

Определите частотный диапазон f и постройте график одностороннего амплитудного спектра P1. Амплитуды находятся не точно в 0,7 и 1, как и ожидалось, из-за добавленного шума. В среднем более длинные сигналы дают лучшие частотные приближения.

f = Fs*(0:(L/2))/L;
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Теперь возьмем преобразование Фурье исходного, некорректированного сигнала и получим точные амплитуды, 0,7 и 1,0.

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)|')

Преобразуйте Гауссов импульс из временного интервала в частотный диапазон.

Задайте параметры сигнала и Гауссов импульс, X.

Fs = 100;           % Sampling frequency
t = -0.5:1/Fs:0.5;  % Time vector 
L = length(t);      % Signal length

X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));

Постройте график импульса во временном интервале.

plot(t,X)
title('Gaussian Pulse in Time Domain')
xlabel('Time (t)')
ylabel('X(t)')

Как использовать fft функция для преобразования сигнала в частотный диапазон, сначала идентифицируйте новую длину входа, которая является степенью 2 от исходной длины сигнала. Это заполнит сигнал X с конечными нулями в порядок улучшить эффективность fft.

n = 2^nextpow2(L);

Преобразуйте Гауссов импульс в частотный диапазон.

Y = fft(X,n);

Задайте частотный диапазон и постройте график уникальных частот.

f = Fs*(0:(n/2))/n;
P = abs(Y/n).^2;

plot(f,P(1:n/2+1)) 
title('Gaussian Pulse in Frequency Domain')
xlabel('Frequency (f)')
ylabel('|P(f)|^2')

Сравните косинусоидные волны в временной интервал и частотный диапазон.

Задайте параметры сигнала с частотой дискретизации 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];

Вывод первых 100 значений из каждой строки X в одной фигуре по порядку и сравните их частоты.

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 аргумент для использования fft вдоль строк X, то есть для каждого сигнала.

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