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

Представление частотного диапазона содержит информацию о величине и фазе сигнала на каждой частоте. Вот почему выход расчета БПФ сложен. Комплексное число,, имеет действительную часть, и мнимую часть, такую, что. Величина вычисляется как, а фаза вычисляется как. Можно использовать функции 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, чтобы восстановить сигнал времени. Флаг 'symmetric' подсказывает 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 кГц непосредственно с выхода БПФ (делая величинам равной нулю) и прослушайте влияние, которое это оказывает на звук аудио файла. Удаление высокой частоты компонентов сигнала называется lowpass фильтрацией.

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

Измерительная Степень

The 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])

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

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

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

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

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

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

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 эффективно удаляет весь ложный peaks частоты, вызванные шумом. Спектральный компонент на 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 - длина наблюдения сигнала. Будут разрешены только спектральные компоненты, разделенные частотой, большей, чем частотное разрешение.

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

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

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

Спектры степени ускорения разомкнутого контура и замкнутого контура показывают, что, когда система управления активна, спектр степени ускорения уменьшается между 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])

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

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

Заключения

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

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

Для получения дополнительной информации об анализе частотного диапазона смотрите Signal Processing Toolbox.

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

Приложение

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