Информация о величине и фазе БПФ

Представление сигнала в частотной области несет информацию о величине и фазе сигнала на каждой частоте. Поэтому вывод вычисления БПФ сложен. Комплексное число,, имеет действительную часть,, и мнимую часть,, что. Величина вычисляется как, а фаза вычисляется как. Можно использовать функции MATLAB. abs и angle чтобы получить соответственно величину и фазу любого комплексного числа.

Используйте звуковой пример, чтобы получить некоторое представление о том, какая информация переносится по величине и фазе сигнала. Для этого загрузите аудиофайл, содержащий 15 секунд акустической гитарной музыки. Частота дискретизации звукового сигнала составляет 44,1 кГц.

Fs = 44100;
y = audioread('guitartune.wav');

Использовать fft для наблюдения за частотным содержанием сигнала.

NFFT = length(y);
Y = fft(y,NFFT);
F = ((0:1/NFFT:1-1/NFFT)*Fs).';

Выход БПФ представляет собой комплексный вектор, содержащий информацию о частотном содержании сигнала. Величина указывает на прочность частотных компонентов относительно других компонентов. Фаза показывает, как все частотные компоненты выравниваются по времени.

Постройте график величины и фазовых составляющих частотного спектра сигнала. Величина удобно наносится на график в логарифмическом масштабе (дБ). Фаза развертывается с помощью unwrap функция, чтобы мы могли видеть непрерывную функцию частоты.

magnitudeY = abs(Y);        % Magnitude of the FFT
phaseY = unwrap(angle(Y));  % Phase of the FFT

helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

Для восстановления временного сигнала можно применить обратное преобразование Фурье к вектору частотной области Y. Флаг «симметричный» сообщает ifft что вы имеете дело с сигналом реального времени, так что он обнуляет небольшие мнимые компоненты, которые появляются при обратном преобразовании из-за численных неточностей в вычислениях. Обратите внимание, что исходный сигнал времени y и восстановленный сигнал y1 практически одинаковы (норма их разности порядка 1e-14). Очень небольшая разница между ними также обусловлена упомянутыми выше числовыми неточностями. Воспроизведение и прослушивание непреобразованного сигнала y1.

y1 = ifft(Y,NFFT,'symmetric');
norm(y-y1)

ans =

   3.9523e-14

hplayer = audioplayer(y1, Fs);
play(hplayer);

Чтобы увидеть эффекты изменения амплитудной характеристики сигнала, удалите частотные составляющие выше 1 кГц непосредственно с выхода БПФ (сделав величины равными нулю) и прислушайтесь к влиянию, которое это оказывает на звук аудиофайла. Удаление высокочастотных составляющих сигнала называется фильтрацией нижних частот.

Ylp = Y;
Ylp(F>=1000 & F<=Fs-1000) = 0;

helperFrequencyAnalysisPlot1(F,abs(Ylp),unwrap(angle(Ylp)),NFFT,...
  'Frequency components above 1 kHz have been zeroed')

Возвращение отфильтрованного сигнала во временную область с помощью ifft.

ylp = ifft(Ylp,'symmetric');

Воспроизвести сигнал. Вы все еще можете услышать мелодию, но это звучит так, как если бы вы закрыли уши (вы фильтруете высокочастотные звуки, когда вы делаете это). Несмотря на то, что гитары производят ноты, которые находятся между 400 и 1 кГц, когда вы играете ноту на струне, струна также вибрирует на кратных базовой частоте. Эти высокочастотные компоненты, называемые гармониками, придают гитаре особый тон. Когда вы их снимаете, вы делаете звук «непрозрачным».

hplayer = audioplayer(ylp, Fs);
play(hplayer);

Фаза сигнала имеет важную информацию о том, когда вовремя появляются ноты песни. Чтобы проиллюстрировать важность фазы звукового сигнала, удалите информацию о фазе полностью, взяв величину каждой частотной составляющей. Обратите внимание, что при этом отклик величины остается неизменным.

% Take the magnitude of each FFT component of the signal
Yzp = abs(Y);
helperFrequencyAnalysisPlot1(F,abs(Yzp),unwrap(angle(Yzp)),NFFT,[],...
  'Phase has been set to zero')

Верните сигнал во временную область и воспроизведите звук. Невозможно распознать исходный звук. Отклик по величине тот же, никакие частотные компоненты на этот раз не были удалены, но порядок нот исчез полностью. Сигнал теперь состоит из группы синусоид, все выровнены в момент времени, равный нулю. В общем случае фазовые искажения, вызванные фильтрацией, могут повредить сигнал до точки его распознавания.

yzp = ifft(Yzp,'symmetric');
hplayer = audioplayer(yzp, Fs);
play(hplayer);

Периодичность поиска сигнала

Представление сигнала в частотной области позволяет наблюдать несколько характеристик сигнала, которые либо нелегко увидеть, либо вообще не видны при взгляде на сигнал во временной области. Например, анализ частотной области становится полезным при поиске циклического поведения сигнала.

Анализ циклического поведения температуры в офисном здании

Рассмотрим набор измерений температуры в офисном здании в зимний сезон. Измерения проводили каждые 30 минут в течение примерно 16,5 недель. Посмотрите на данные временной области с временной осью, масштабированной до недель. Может ли быть какое-либо периодическое поведение в отношении этих данных?

load officetemp.mat
Fs = 1/(60*30);   % Sample rate is 1 sample every 30 minutes
t = (0:length(temp)-1)/Fs;

helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp,...
  'Time in weeks','Temperature (Fahrenheit)')

Практически невозможно узнать, существует ли какое-либо циклическое поведение на рабочих температурах, глядя на сигнал временной области. Однако циклическое поведение температуры становится очевидным, если посмотреть на ее представление в частотной области.

Получение представления сигнала в частотной области. Если построить график величины выхода БПФ с частотной осью, масштабированной до циклов/недель, можно увидеть, что существуют две спектральные линии, которые явно больше любой другой частотной составляющей. Одна спектральная линия лежит на 1 цикле/неделе, другая - на 7 циклах/неделе. Это имеет смысл, учитывая, что данные поступают из контролируемого температурой здания по 7-дневному календарю. Первая спектральная линия указывает на то, что температура зданий следует недельному циклу с более низкими температурами в выходные дни и более высокими температурами в течение недели. Вторая линия указывает на то, что существует также дневной цикл с более низкими температурами ночью и более высокими температурами днем.

NFFT = length(temp);              % Number of FFT points
F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector

TEMP = fft(temp,NFFT);
TEMP(1) = 0; % remove the DC component for better visualization

helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)),...
  'Frequency (cycles/week)','Magnitude',[],[],[0 10])

Измерительная мощность

periodogram функция вычисляет БПФ сигнала и нормализует выходной сигнал для получения спектральной плотности мощности, PSD или спектра мощности, по которому можно измерить мощность. PSD описывает, как мощность временного сигнала распределяется с частотой, она имеет единицы мощности/Гц. Вы вычисляете спектр мощности, интегрируя каждую точку PSD на частотном интервале, на котором эта точка определена (т.е. на разрешающей полосе полосы PSD). Единицами спектра мощности являются ватты. Значения мощности можно считывать непосредственно из спектра мощности без необходимости интегрироваться через интервал. Обратите внимание, что PSD и спектр мощности являются реальными, поэтому они не содержат никакой фазовой информации.

Измерение гармоник на выходе нелинейного усилителя мощности

Загрузить данные, измеренные на выходе усилителя мощности, имеющего искажение третьего порядка формы, где - выходное напряжение и - входное напряжение. Данные собирали с частотой выборки 3,6 кГц. Вход состоит из синусоиды 60 Гц с единичной амплитудой. Из-за характера нелинейных искажений следует ожидать, что выходной сигнал усилителя будет содержать постоянную составляющую, составляющую 60 Гц и вторую и третью гармоники на частоте 120 и 180 Гц.

Загрузка 3600 выборок выходного сигнала усилителя, вычисление спектра мощности и построение графика результата в логарифмическом масштабе (децибелы-ватты или дБВт).

load ampoutput1.mat
Fs = 3600;
NFFT = length(y);

% Power spectrum is computed when you pass a 'power' flag input
[P,F] = periodogram(y,[],NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

График спектра мощности показывает три из четырех ожидаемых пиков при DC, 60 и 120 Гц. Он также показывает еще несколько ложных пиков, которые должны быть вызваны шумом в сигнале. Заметим, что гармоника 180 Гц полностью погружена в шум.

Измерьте мощность видимых ожидаемых пиков:

PdBW = 10*log10(P);
power_at_DC_dBW = PdBW(F==0)   % dBW

[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,'minpeakheight',-11);
peakFreqs_Hz = F(peakFreqIdx)
peakPowers_dBW

power_at_DC_dBW =

   -7.8873


peakFreqs_Hz =

    60
   120


peakPowers_dBW =

   -0.3175
  -10.2547

Улучшение измерений мощности для шумных сигналов

Как видно на графике выше, периодограмма показывает несколько пиков частоты, которые не связаны с интересующим сигналом. Спектр выглядит очень шумным. Причина этого в том, что вы проанализировали только одну короткую реализацию шумного сигнала. Повторение эксперимента несколько раз и усреднение удалят ложные спектральные пики и дадут более точные измерения мощности. Вы можете достичь этого усреднения с помощью pwelch функция. Эта функция будет принимать большой вектор данных, разбивать его на меньшие сегменты заданной длины, вычислять столько периодограмм, сколько есть сегментов, и усреднять их. По мере увеличения числа доступных сегментов pwelch функция даст более плавный спектр мощности (меньшая дисперсия) с значениями мощности, более близкими к ожидаемым значениям.

Загрузите большее наблюдение, состоящее из 500e3 точек выхода усилителя. Количество точек, используемых для выполнения FFT, должно быть равно 3600, так что floor (500e3/3600) = 138 FFT усредняются для получения спектра мощности.

load ampoutput2.mat
SegmentLength = NFFT;

% Power spectrum is computed when you pass a 'power' flag input
[P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
  'Power spectrum (dBW)',[],[],[-0.5 200])

Как видно на сюжете, pwelch эффективно удаляет все ложные пики частоты, вызванные шумом. Теперь видна спектральная составляющая при 180 Гц, которая была погружена в шум. Усреднение удаляет дисперсию из спектра, и это эффективно дает более точные измерения мощности.

Измерение общей средней мощности и мощности в диапазоне частот

Измерение общей средней мощности сигнала временной области является простой и общей задачей. Для выходного сигнала усилителя, y, общая средняя мощность вычисляется во временной области как:

pwr = sum(y.^2)/length(y) % in watts

pwr =

    8.1697

В частотной области общая средняя мощность вычисляется как сумма мощности всех частотных составляющих сигнала. Значение pwr1 состоит из суммы всех частотных составляющих, имеющихся в спектре мощности сигнала. Значение согласуется со значением pwr, вычисленным выше с использованием сигнала временной области:

pwr1 = sum(P) % in watts

pwr1 =

    8.1698

Но что, если вы хотите измерить общую мощность, доступную в диапазоне частот? Вы можете использовать bandpower функция для вычисления мощности по любой требуемой полосе частот. Вы можете передать сигнал временной области непосредственно в качестве входа в эту функцию, чтобы получить мощность по указанному диапазону. В этом случае функция будет оценивать спектр мощности методом периодограммы.

Вычислите мощность в диапазоне от 50 Гц до 70 Гц. Результат будет включать в себя мощность 60 Гц плюс мощность шума в интересующей полосе:

pwr_band = bandpower(y,Fs,[50 70]);
pwr_band_dBW = 10*log10(pwr_band) % dBW

pwr_band_dBW =

    0.0341

Если требуется управлять вычислением спектра мощности, используемого для измерения мощности в полосе, можно передать вектор PSD в bandpower функция. Например, можно использовать pwelch функция, как и прежде, для вычисления PSD и обеспечения усреднения шумовых эффектов:

% Power spectral density is computed when you specify the 'psd' option
[PSD,F]  = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'psd');
pwr_band1 = bandpower(PSD,F,[50 70],'psd');
pwr_band_dBW1 = 10*log10(pwr_band1) % dBW

pwr_band_dBW1 =

    0.0798

Поиск спектральных компонентов

Сигнал может состоять из одной или более частотных составляющих. Возможность наблюдения за всеми спектральными компонентами зависит от частотного разрешения анализа. Частотное разрешение или разрешающая полоса спектра мощности определяется как R = Fs/N, где N - длина наблюдения сигнала. Разрешаются только спектральные компоненты, разделенные частотой, превышающей разрешение частоты.

Анализ системы управления вибрациями при землетрясениях в здании

Системы управления Active Mass Driver (AMD) используются для снижения вибрации в здании при землетрясении. На верхнем этаже здания размещается активный массовый драйвер, и на основании измерений смещения и ускорения этажей здания система управления посылает водителю сигналы, чтобы масса двигалась для ослабления возмущений грунта. Измерения ускорения регистрировались на первом этаже трехэтажного испытательного сооружения в условиях землетрясения. Измерения проводились без активной системы управления массовым приводом (состояние разомкнутого контура) и с активной системой управления (состояние замкнутого контура).

Загрузите данные ускорения и вычислите спектр мощности для ускорения первого этажа. Длина векторов данных составляет 10e3, а частота дискретизации - 1 кГц. Использовать pwelch с сегментами длиной 64 точки данных для получения floor (10e3/64) = 156 средних значений БПФ и пропускной способности разрешения Fs/64 = 15,625 Гц. Как было показано ранее, усреднение уменьшает шумовые эффекты и дает более точные измерения мощности. Используйте 512 точек БПФ. Использование NFFT > N эффективно интерполирует частотные точки, формируя более подробный график спектра (это достигается добавлением нулей NFFT-N в конце временного сигнала и взятием NFFT-точки FFT нулевого дополненного вектора).

Спектры мощности ускорения с разомкнутым контуром и с замкнутым контуром показывают, что когда система управления активна, спектр мощности ускорения уменьшается между 4 и 11 дБ. Максимальное ослабление происходит примерно при 23,44 кГц. Уменьшение на 11 дБ означает, что мощность вибрации уменьшается в 12,6 раза. Суммарная мощность снижается с 0,1670 до 0,059 Вт, коэффициент 2,83.

load quakevibration.mat

Fs = 1e3;                 % sample rate
NFFT = 512;               % number of FFT points
segmentLength = 64;       % segment length

% open loop acceleration power spectrum
[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');

% closed loop acceleration power spectrum
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
  'Frequency in Hz','Acceleration Power Spectrum in dB',...
  'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100])

Вы анализируете данные вибрации и знаете, что вибрации имеют циклическое поведение. Тогда каким образом графики спектра, показанные выше, не содержат каких-либо острых спектральных линий, типичных для циклического поведения? Может быть, вы пропускаете эти линии, потому что они не разрешимы с разрешением, полученным с длиной сегмента 64 точки? Увеличьте частотное разрешение, чтобы узнать, существуют ли спектральные линии, которые ранее не разрешались. Для этого увеличьте длину сегмента данных, используемого в pwelch функция до 512 точек. Это дает новое разрешение Fs/512 = 1,9531 Гц. В этом случае число средних значений БПФ уменьшается до пола (10e3/512) = 19. Очевидно, что существует компромисс между количеством средних значений и разрешением частоты при использованииpwelch. Количество точек БПФ должно быть равно 512.

NFFT = 512;            % number of FFT points
segmentLength = 512;   % segment length

[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');
P1_CL     = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');

helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
  'Frequency in Hz','Acceleration Power Spectrum in dB',...
  'Resolution bandwidth = 1.95 Hz',{'Open loop', 'Closed loop'},[0 100])

Обратите внимание, как увеличение разрешения частоты позволяет наблюдать три пика на спектре разомкнутого контура и два на спектре замкнутого контура. Эти пики раньше не разрешались. Расстояние между пиками в спектре разомкнутого контура составляет около 11 Гц, что меньше, чем разрешение по частоте, полученное с сегментами длиной 64, но больше, чем разрешение, полученное с сегментами длиной 512. Теперь видно циклическое поведение вибраций. Частота основных колебаний составляет 5,86 Гц, и пики частоты в равноудаленном состоянии позволяют предположить, что они являются гармонически связанными. Хотя уже было замечено, что система управления уменьшает общую мощность колебаний, спектры с более высоким разрешением показывают, что другим эффектом системы управления является вырезание гармонической составляющей при 17,58 Гц. Таким образом, система управления не только уменьшает вибрацию, но и приближает ее к синусоиде.

Важно отметить, что частотное разрешение определяется количеством точек сигнала, а не количеством точек БПФ. Увеличение числа точек БПФ интерполирует частотные данные для получения более подробной информации о спектре, но это не улучшает разрешение.

Заключения

В этом примере вы научились выполнять анализ частотной области сигнала с помощью fft, ifft, periodogram, pwelch, и bandpower функции. Вы поняли сложный характер БПФ и информацию, содержащуюся в величине и фазе частотного спектра. Вы видели преимущества использования данных частотной области при анализе периодичности сигнала. Вы узнали, как вычислить полную мощность или мощность по определенной полосе частот шумного сигнала. Вы поняли, как увеличение частотного разрешения спектра позволяет наблюдать близко расположенные частотные компоненты, и вы узнали о компромиссе между частотным разрешением и спектральным усреднением.

Дальнейшее чтение

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

Справка: Дж. Г. Проакис и Д. Г. Манолакис, "Цифровая обработка сигналов. принципы, алгоритмы и приложения, "Прентис Холл, 1996.

Приложение

В этом примере используются следующие вспомогательные функции.