Can we differentiate between observed sequences of IID and uncorrelated noise? Given a realization $\{z_{i}\}_{i}$ of some noise process $\{Z_{i}\}_{i}$ that is either a White-Noise process or an IID process, i.e. both are uncorrelated however the former also is independent, can we determine which one it is?

It is not clear what exactly the OP understands or means by the various words in the phrase

either a White-Noise process or an IID process, i.e. both are uncorrelated however the former also is independent,

For example, the initialism IID usually stands for "independent identically distributed", and so "IID process" usually means a sequence of independent identically distributed random variables. If readers agree that this is the accepted meaning, then the OP's statement that the random variables constituting the process are uncorrelated is perfectly true, but the very next words "former also (emphasis added) is independent" suggests that the OP thinks that the random variables constituting the IID process are not independent, uncorrelated though they happen to be. Of course, those who don't accept the generally accepted meaning of IID will have their own version of what is the truth and what is a BIG LIE.
Similarly, it is not clear what the OP understands by "White-Noise process", except that the random variables constituting the process are independent random variables (and hence uncorrelated random variables). Are there any other requirements that might be reasonable to impose on the random variables constituting the White-Noise process? Zero mean, maybe? Symmetric distribution about $0$? Identical Distribution? or is something like $X[n] \sim \mathcal N(0,1)$ for even $n$ and $X[n] \sim \mathcal N(0,2)$  for odd $n$ acceptable? And if so, what is the justification for this belief about the properties of a White Noise process?
There is a long-standing tradition on dsp.SE to say that (discrete-time) white noise is nothing more than a sequence of uncorrelated random variables and that independence is not required, and should not be insisted upon. I disagree with this school of thought because (i) correlation is a pairwise relation whereas independence is a more general property involving more than two random variables and (ii) lack of correlation only tells us that linear estimation/prediction is not useful but says nothing about more general nonlinear estimators/predictors.  A canonical example is of $X
\sim \mathcal N(0,1)$ and $Y = X^2$ which are uncorrelated random variables and yet $Y$ is a (nonlinear) function of $X$. As another example. consider the following from @MarcusMuller's answer to this question:

Dependent, but uncorrelated (in the pair-wise sense)
Also, $Y_i$ be a discrete time process. Define:
$$Y_i \begin{cases}
>\sim U\{-1,+1\} &\text{randomly, independently}&  \text{ for }i \ne 3N,\\
>=Y_{i-1}\cdot Y_{i-2} &\text{(fixed!)} & \text{ for }i=3N
>\end{cases}.$$

Indeed, $Y_1, Y_2,$ and $Y_3 = Y_1Y_2$ uncorrelated random variables, but in fact, $Y_1, Y_2, Y_3$ are pairwise independent random variables too.  However, $Y_1, Y_2, Y_3$ are not mutually independent random variables.


check power spectrum my question is quite straight-forward to be grasped and handled.

I have an 1D array, F_11 representing a target power spectrum. It is sketched in the following figure

Then, I start generating a time series u which embodies the properties of F_11, as follows

nn = complex(normrnd(0,1,[1,1000]),normrnd(0,1,[1,1000]));
u = abs(fftn(fftshift(F_11.*nn)));

Finally, just to check that u does actually reproduce a spectrum as F_11, I try going backwards as follows

ps =  abs(ifftn(ifftshift(u)));

But, when plotting ps, it does not exactly match F_11

In the high frequency area, ps drifts upwards, which is quite weird.

Would you mind to sketch the flaw in my procedure?

Your u variable is just the magnitude of the FFT, not the actual complex numbers that represent the FFT. Try removing the abs() from u and then doing the IFFT.


What does Fast and Slow Time Weighting signify for real-time audio signal processing? I have a streaming audio signal that I am recording in real time. To make the explanation easy, say the signal that I receive is calibrated and has a sampling frequency, $F_s = 44100\textrm{ Hz}$.

I am then finding the Sound Pressure Level of this signal. For this purpose, I get a frame of audio data of t seconds.

Then, RMS Sound Pressure Level (dB) is given by:
$$
\rm Amplitude_{RMS} = \sqrt{mean\left(frame ^2\right)}
$$

This gives the instantaneous RMS of a $t$-second frame of audio. Then, to find a $\textrm{ dB}$ level,
$$
L_p = 20\log_{10}\left(\rm Amplitude_{RMS}\right)
$$

Things get tricky here with the definitions of 'Fast' or 'Slow' sound pressure level. As this post points out, for time-weighted exponential average, I created a low-pass filter with a real pole at $1/\tau$ (where $\tau= 125\textrm{ ms}$ for 'fast' or $1\textrm{ s}$ for 'slow')

τ = 0.125; %fast
[b a] = bilinear(1, [1 1/τ], Fs);
frame_filtered = filter(b, a, frame.^2);
L_fast = 20*log10(sqrt(mean(frame_filtered)));

The IEC 61672 specifies the time weightings to be Fast and Slow for determining the speed at which the instrument responds to changing noise levels.

How does this work for digital audio?
Does it mean that the analysis needs to be performed on a fixed frame of data?
For this particular example, this frame would be $t\cdot F_s$ samples?
Slow: 1*44100 = 441000 samples
Fast: 0.125*44100 ~= 5513 samples
If so, is it true that I would need to know in advance whether I need fast or slow levels?

By "time weighting" in acoustics we are applying an RC (resistor–capacitor) circuit to a time signal.

See: RC circuit on Wikipedia

The time constant for FAST is 0,125 s and 1 s for SLOW.

You can also apply this digitally on your time signal by using my code below - it is basically a exponential moving average filter. I think you understand the steps by just looking at the comments in the code below.

  public class FilterTimeWeighting {

private int nBands = 29; //number of bands
private float initValue;

FilterTimeWeighting() {
    initValue=0;
}

/**
 * Apply the LAmax RC-filter on the signal. Valid for fsamp=44100 Hz.
 * Returns the signal in Volts after applying filter.
 *
 * //y(n)=a*y(n-1)+(1-a)*x(n)
 *
 *
 *
 * @param signalIn
 * @return
 */
public float[] applyLmaxFilter(float[] signal) {

    int n_samples = signal.length;

    float[] x_signal = new float[n_samples];
    float[] y_signal = new float[n_samples];
    float[] ret_signal = new float[n_samples];

    //Square the signal
    for (int i_sample = 0; i_sample < n_samples; i_sample++) {
        x_signal[i_sample] = signal[i_sample] * signal[i_sample];
    }

    //Exponential moving average parameters
    float tau=0.125f;
    float fs=44100f;
    float alpha=tau/(1f/fs+tau);

    //for the first sample
    y_signal[0] = initValue*alpha+(1-alpha)*x_signal[0];

    //for the rest of the samples
    for (int i_sample = 1; i_sample < n_samples; i_sample++) {
        y_signal[i_sample] = y_signal[i_sample-1]*alpha+(1-alpha)*x_signal[i_sample];
    }

    //remember last sample for next time step;
    initValue=y_signal[n_samples-1];

    //Square root
    for (int i_sample = 0; i_sample < n_samples; i_sample++) {
        ret_signal[i_sample] = (float) Math.sqrt(y_signal[i_sample]);
    }

    return ret_signal;
}

  }


Extracting frequencies from FFT I performed 512 point FFT on a signal. I got another set of  512 Numbers. I understand that those numbers represent amplitude of the various sine and cosine waves having different frequencies.

If my understanding is correct can somebody tell me how to know the frequencies of those sine and cosine waves from the knowledge of those 512-numbers (i.e amplitudes) ?

Assuming your 512 samples of the signal are taken at a sampling freqeuncy $f_s$, then the resulting 512 FFT coefficients correspond to frequencies:

0, $f_s/512$, $2 f_s/512$, $\ldots$, $511 f_s/512$

Since you are dealing with discrete-time signals, Fourier transforms are periodic, and FFT is no exception.

Therefore you can think of your last coefficient as corresponding also to frequency $511 f_s/512 = (511-512) f_s/512 = -1 f_s/512$.

The same applies to the second to last coefficient, and so on. This is the mirroring commented by Daniel Hicks.

Also, if you are transforming a real signal, then all your information is contained in the first 256 FFT coefficients. The rest are simply complex conjugates of the first coefficients.


AutoCorrelation using FFT of power spectrum I have a chunk of data (samples) that need to do autocorrelation using FFT method. All data are real and lies between -1 and +1. What I have done follows:

Zero pad the window with the length equals the samples' length.
FFT of result from 1
Result from 2 has Real and Imaginary parts. I used them to find Power Spectral Density by taking Re*Re + Im*Im.
Take the inverse Fast Fourier Transform.

Now here is my question: The result from part 4 contains Real and Imaginary parts. Is the magnitude of it equals the autocorrelation result? If so, since my input contains positive and negative values, the magnitude of Result from part 4 means the autocorrelation result cannot be negative, which doesn't make sense.

Or should I take only the Real part of result from Part 4 as the autocorrelation result? My input is real and so the autocorrelation shouldn't have any imaginary part, right?

In principle you are doing the right thing. Step 4 should produce a real result unless there is a coding error. Sometimes you add up some residual imaginary part due to numerical noise but, if any, that should be very small.

Here is an example

%% random signal of length n
n = 128;
x = rand(n,1);
% zero pad
x = [x; zeros(n,1)];
% fft
fx = fft(x);
% mag sqaured
fa = fx.*conj(fx);
% inverse FFT
a = ifft(fa);
% plot
plot(-n+(1:2*n),circshift(a,-n));
answer = {'no','yes'};
fprintf('Is a Real?: %s \n',answer{isreal(a)+1});


Octaves vs Decades Fractions At this USGS very good site about rms noise in sensors, a phrase regarding third octave bands (the original topic) appeared.

  The figure displays the minimum vertical seismic noise directly as rms amplitudes in a bandwidth of one-sixth decade. By coincidence, these amplitudes may also be interpreted as average peak amplitudes in a bandwidth of one-third octave (a standard bandwidth in acoustics)

One sixth decade factor is $10^{1/6}=1.467799$; the 10th factor $10/6=1.666666$, and one third octave factor $2^{1/3}=1.259921$, none of them being comparable.

It happens actually to be comparable with one tenth decade factor $10^{1/10}=1.258925$.

Note the rms and peak distictions... Why the author state as coincident both log factors?

In addition, the same paragraph continues with:

  An example: the minimum vertical ground noise between the periods of 10 and 20 sec is at -180 dB relative to  ${\rm 1 m/s^2}$, thus  ${\rm 10^{-180/20}\;m/s^2 = 1\;nm/s^2}$ average peak in one-third octave.

Which is ok, but then:

  Since the bandwidth considered is a full octave, the total average peak amplitude in this band is  ${\rm\sqrt{3}\;nm/s^2}$.

In here i honestly dont know why a $\sqrt 3$ factor appeared...

I think they made a mistake. $2^{1/3} \approx 1.259921$ is not close enough to $10^{1/6} \approx 1.467799$ to say that there is any coincidence. Even the half-octave factor $2^{1/2} \approx 1.414214$ is closer.


Questions on Need of filtering before decimation? This is a question related to the need for filtering before decimation. I understand that there are many questions likewise in this forum. I have gone through these and still little confused (or) may be overthinking about it.

I have a test sequence(1x1000)sampled at 176.4 KHz(example) and I need to decimate it to let's say 44.1KHz. (the maximum frequency component in the sequence is 20KHz, I am taking examples from audio freqs).
The way I understand decimation is, any sequence with sample rate Fs is decimated to Fsd(=Fs/D) ... then there is a higher chance that it violates the Nyquist sampling theorem and information is lost in the original sequence(if we don't think about filtering and just decimate it).
In my example, if I decimate (by factor 4) 176.4Khz to 44.1Khz, it still satisfies the Nyquist sampling, so there will be no concept of aliasing since Fsd>2*20Khz ------ **In this case, do I need to apply LPF before Decimation. If Yes why ?*
If I want to decimate by factor 16, i.e. 176.4KHz to 11.025Khz, it doesn't satisfy Nyquist sampling (as my new Fsd<2*20KHz).. aliasing occurs.So, we (undesirably) remove the higher frequency components in my test sequence by applying LPF (with cutoff freq ≤ 5.515 (F_coutoff ≤ Desired_sample_rate/2) )

LPF is applied such that there are no frequency components higher than 5.515KHz
Is this how I need to choose my cut-off ??

If you think for a moment why you filter, it is to remove frequencies that would fold to the remaining bandwidth when decimating.

So if these frequencies that would fold do not exist in the audio then it will have same bandwidth whether you filter it or not. So in the first case, there is really no need to apply LPF first (unless there is noise that would fold and cause problems).

In the second case, it needs to be filtered, and you calculated the cut-off correctly. However, in practice it won't be a brick-wall filter, so depending on needed filter performance and allowable transition bandwith, you might want to select for example cut-off of 4.5 kHz and allow the filter to have about 1kHz transition band to have some reasonable attenuation level at 5.515 kHz.


Relation between time domain, DTFT domain and frequency domain Problem

The sampling frequency of a continuous-time signal is $S$ kHz, what does $\frac{\pi}{4}$ radians/sample in DTFT domain represent in Hz in frequency domain? Prove the relationship.

Doubts

I am confused with DTFT domain and frequency domain and I do not how to put all these domains together.

If this is a homework problem, shield your eyes.

Otherwise, it is a simple units issue.

S has units of 1k samples per second.

Your frequency is expressed in radians per sample.

There are $ 2 \pi $ radians per cycle.

So the puzzle now is to find the relationship between $\frac{samples}{second}$, $\frac{radians}{sample}$, $\frac{radians}{cycle}$, and your goal of $\frac{cycles}{second}$,


Estimate input from measuring only RC filtered input I am interested in knowing the step amplitude on a signal which can only be measured through a hardware RC filter. The measured signal follows this formula:

$$
V(t) = V_0 \left(1 - e^{-\frac{t}{RC}}\right)
$$

Is it possible to estimate the step amplitude, $V_0$, given that we can perform two measurements on the almost linear part of the rising edge (with a known delay between them)?

That's certainly possible, assuming you know the values of the resistance and capacitance $R$ and $C$. You can expand the signal $V(t)$ in a Taylor series about $t=0$ to arrive at an approximate linear model. First, remember the Taylor series of the exponential function:

$$
e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!} = 1 + x + \frac{x^2}{2} + \ldots
$$

If you're assuming a linear model near the step input time of $t=0$, then truncate the sum to arrive at:

$$
e^x \approx 1 + x
$$

Substitute this into the equation for $V(t)$ to yield:

$$
V(t) \approx V_0\left(1 - \left(1 - \frac{t}{RC}\right)\right)
$$

$$
V(t) \approx V_0\frac{t}{RC}
$$

Which is a pretty simple result. If you know the absolute time of the step input, then you theoretically don't even need multiple measurements; you could estimate $V_0$ directly using:

$$
V_0 \approx \frac{V(T)RC}{T}
$$

where $T$ is the time at which you measure the signal. If you don't know the exact time at which the step is input, then you can take two measurements, as you suggested, at times $T_1$ and $T_2$, then form the estimate as follows:

$$
V(T_2) - V(T_1) \approx V_0\frac{T_2}{RC} - V_0\frac{T_1}{RC}
$$

$$
V_0 \approx RC \frac{V(T_2) - V(T_1)}{T_2 - T_1}
$$


Hilbert transform pair proof I am looking for the proof that the Hilbert transform of $\displaystyle\frac{\sin(at)}{at}$ is given by

$$\frac{\sin^2(at/2)}{at/2}.$$

How do we prove this? This is a $\operatorname{sinc}(at)$ function. I could not find a proof of this in any textbook so far. Preferably I want to arrive at this result using the analytic signal whose imaginary part would be the Hilbert transform:

$$f_a(t) = 2 \int^\infty_0 FT \Big[ \frac{\sin(at)}{at} \Big] \exp(j2 \pi \nu t) \ d\nu$$

I tried this in a few different ways, but could not get to the answer.

I agree that one of the easiest ways to compute the Hilbert transform in this case is to use the analytic signal. This is most easily obtained via the Fourier transform. Note that the Fourier transform of

$$f(t)=\frac{\sin(at)}{at}\tag{1}$$

is

$$F(j\omega)=\frac{\pi}{a}\text{rect}\left(\frac{\omega}{2a}\right)\tag{2}$$

i.e., a rectangular function with value $\pi/a$ in the interval $[-a,a]$, and zero otherwise. The spectrum of the corresponding analytic signal equals the part of $F(j\omega)$ at positive frequencies, scaled by a factor of $2$:

$$F_a(j\omega)=\frac{2\pi}{a}\text{rect}\left(\frac{\omega}{a}-\frac12\right)\tag{3}$$

The inverse Fourier transform of $(3)$ gives the analytic signal:

$$f_a(t)=\frac{1}{2\pi}\frac{2\pi}{a}\int_0^{a}e^{j\omega t}d\omega=\frac{1}{jat}\left(e^{jat}-1\right)=e^{jat/2}\frac{\sin(at/2)}{at/2}\tag{4}$$

The Hilbert transform of $f(t)$ is the imaginary part of $(4)$:

$$\hat{f}(t)=\text{Im}\{f_a(t)\}=\frac{\sin^2(at/2)}{at/2}\tag{5}$$

as stated in your question.


How Can I Remove Flickering Background in a Video? I have a black and white video taken with a high speed camera and I notice a background flickering coming from my LED light.  The images alternate from dim to bright and it makes the videos quite disorienting.

Besides getting a different light, is there any post-proccessing that can be done to remove this effect?  Possibly something from the image processing toolbox in matlab, or something in imageJ would be helpful as well.

I have already acquired quite a few videos and would hate to have to repeat everything all over again so any help would be greatly appreciated.

thanks!

For those classic Video Processing operations there is nothing better than the Plug In's of AviSynth.

Specifically for De Flickering look at:

LMFlicker.
ReduceFlicker.
DeFlicker.

The source code of each is available so you can get inspiration from them.

All of them basically works on smoothing the Luma Data along the temporal dimension.


Power Spectral Density of discrete White Gaussian Noise defined with variance and sampling interval From Creighton and Anderson Chapter 7 on Gravitational Wave data analysis:

Please explain where the delta t and the limit came from in this derivation.

Marcus is correct, so this means something else is wrong. The left equality simply shows that, for continuous time, the unilateral PSD, $S_x(f)$, is two times the Fourier transform of the autocorrelation function, $R_x(\tau)$. See, e.g., A.B. Carlson, Introduction to Communication Systems, 2nd Ed., McGraw-Hill, NY ©1975, Chapter 2, equation 18a. Or see pretty much any other book on Fourier transforms.

In discrete time, with consecutive samples spaced by $\Delta t$, the unilateral PSD for white noise is $2\sigma^2\Delta t$. Small values of $\Delta t$ are usually desirable, since the sampling frequency is the reciprocal of $\Delta t$ and the Nyquist frequency is half the sampling frequency. But $\Delta t$ cannot go to zero, as Marcus commented.

What happens if you just delete the limit stuff, on the right hand side of equation 7.22, and see what happens in the textbook's subsequent matched filter treatment? For a matched filter, the maximum signal-to-noise ratio is simply twice the pulse energy divided by the unilateral (white noise) PSD. (Same reference as above, Chapter 4, equation 12a.) So it looks like $2\sigma^2\Delta t$ will be used in your textbook. Hope this helps a bit.


How many moments are needed for a discrete distributions I am doing a project on imaging and I have a more general question about imaging.

In general we use the moments of the distribution of the object imaged, because the capture information about the image. If our image is continuous then we will need in general infinite number of moments to describe it correctly. (Because if we consider the image as a function and we take the Taylor expansion in moments we will need all of them).

My question is, if we have a discrete distribution, 2 or 5 points that we want to image, how many moments do we need? Do we still need infinite of them or just 2 or 5?

Thank you!

You can do $E\{ x^n \}$ as much as you want. How much is necessary is up to you but the unique characterization of a distribution is not the same thing as a cutoff point of a moment expansion.  Some distributions can be very compact with respect to their moments but others are not.  My experience is that moment (taylor like) expansions can converge very slowly.

For something like a binomial distribution, all I need to know to represent the pdf is $p=$ success and know that the trials are independent.   The probability density for successes given a number of trials flollows from the mathematical formula.


Possible to implement simple IIR filter with greater than 6 db/octave slope? I am using the following IIR code as a low-pass filter in a software project, and it works great. However, I would like the slope of the filter to be much steeper as the one here seems to be very gradual (6 db/octave, from what I understand):

#define CUTOFF_FREQ 500.0
#define SAMPLE_RATE 48000.0

float RC = 1.0/(CUTOFF_FREQ*2*M_PI);
float dt = 1.0/SAMPLE_RATE;
float alpha = dt/(RC+dt);

float output_sample,prev_sample,cur_sample;

/* Run the filter on the two samples (normally done inside a loop) */
output_sample = prev_sample + (alpha*(cur_sample - prev_sample));

I would like perhaps a 12 db/octave slope or even 18 db/octave, similar to the low-pass filters in popular audio programs. Is it possible to modify this code for this purpose? Or should I be looking into alternate implementations altogether?

You need a filter with a higher filter order, which means that your filter memory must be longer. The current output sample $y[n]$ is computed as

$$y[n]=a_1y[n-1]+a_2y[n-2]+\ldots + a_Ny[n-N]+b_0x[n]+b_1x[n-1]+\ldots + b_Nx[n-N]$$

where $x[n]$ is the current input sample, $N$ is the filter order, and $a_i$ and $b_i$ are the filter coefficients. You can design the filter (i.e., compute the coefficients) using some tool like Matlab or Octave.


Modulating ofdm in audio spectrum Trying to modulate audio with OFDM data.

My fs = 44.1khz, i'm using 64 sub-carriers, but i'm only using 5 sub-carriers for data, so my bandwidth is about 3.445khz. I'm not sure how to properly write the all real data to a wav file and then reconstruct the complex signal when I read it back. Can someone point me where I'm messing up?

Matlab code:

 close all
clear all

%fft length
N=64;
NN=1024;

%carrier frequency
fc = 15000;

%create the data to send
qpsk1 = (floor(2*rand(1,2))-.5)/.5 + 1j*(floor(2*rand(1,2))-.5)/.5;
qpsk2 = (floor(2*rand(1,3))-.5)/.5 + 1j*(floor(2*rand(1,3))-.5)/.5;

inputiFFT = [zeros(1,29) qpsk1, 0, qpsk2 zeros(1,29)];
outputiFFT = ifft(inputiFFT,N);

%spin data to carrier
y3 = outputiFFT.*exp(1j*2*pi*15000/44100*(0:length(outputiFFT)-1));
%despin to baseband for sanity check
y4 = y3.*exp(-1j*2*pi*15000/44100*(0:length(y3)-1));

figure
subplot(211)
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(outputiFFT,NN)))))
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))))
hold off
title('Data - blue, Data at carrier - r')
subplot(212)
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y4,NN)))))
title('sanity check down mix')

%upconvert to 15khz and write to wav (all 'real')
towav = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)) + imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

%write to wav then read from wave
wavwrite(towav, 44100, 'ofdmTest.wav');
[y_wav,FS] = wavread('ofdmTest.wav');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(towav,NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','r')
hold off
title('Theoretical data up converted to carrier -r , data written to wav - b')

%convert to I-Q
yRx = (y_wav.*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)).'+ 1j*y_wav.*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1)).');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(yRx.',NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y_wav,NN)))),'color','r')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','g')

hold off
title('yRx -b, data written to wav -r, theoretical - g')

Typically you use Orthogonal Frequency Division Multiplexing (OFDM) with communications. In this context having complex data in the time domain makes sense because the real and imaginary parts of the number represent the inphase and quadrature parts of the signal.

In an audio signal which is what a wav file is used for is only real. As sound is a naturally occurring signal.

I am not sure why you'd want to store the inverse of an FFT in a wav file (listening to it would sound awful). As you've discovered there's no way of guaranteeing that this data is real. You seem to get around this by taking the real part of the output and only writing that to your file.

However, if I assume you need this data stored in a wav file you could store the real and imaginary part in the wav file by interleaving them. You could write something like:

towav = zeros(length(outputiFFT)*2);
towav(1:2:end-1) = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1));
towav(2:2:end) = imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

in your code above.

Bear in mind this "halves” your sample rate when you write the 'towav' variable to your wav file. Although I am assuming you don't need to listen to the file, so this is unlikely to matter.

EDIT: as the comment below points out if you weren’t generating random QPSK signals and you ensured that the input sequence in the frequency domain was conjugate symmetric you would have a real sequence in the time domain and you would not need to store the imaginary part.


PSD scaling due to the windowing I am struggling with the following problem for a while: If I have signal of duration N and I apply window (assume rectangular window) of size M < N how will this influence PSD of that signal after doing fft?

PSD of signal should be scaled (I assume) but by which amount and what is the mathematical explanation behind it?

Thank you for your answers in advance.

The shortest answer to that question is to divide your PSD by the squared sum of window values $w_i$:

$$PSD(n) = \dfrac{2|X(n)|^2}{S^2} $$

where:

$n=0\ldots N/2 $

$S = \sum_{i=0}^{N-1}w_i$

You can read more on that topic in the following article:

  Heinzel G., et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows.


What is the cut-off frequency of a moving average filter? I need to design a moving average filter that has a cut-off frequency of 7.8 Hz. I have used moving average filters before, but as far as I'm aware, the only parameter that can be fed in is the number of points to be averaged... How can this relate to a cut-off frequency?

The inverse of 7.8 Hz is ~130 ms, and I'm working with data that are sampled at 1000 Hz. Does this imply that I ought to be using a moving average filter window size of 130 samples, or is there something else that I'm missing here?

The moving average filter (sometimes known colloquially as a boxcar filter) has a rectangular impulse response:

$$
h[n] = \frac{1}{N}\sum_{k=0}^{N-1} \delta[n-k]
$$

Or, stated differently:

$$
h[n] = \begin{cases}
\frac{1}{N}, && 0 \le n < N \\
0, && \text{otherwise}
\end{cases}
$$

Remembering that a discrete-time system's frequency response is equal to the discrete-time Fourier transform of its impulse response, we can calculate it as follows:

$$
\begin{align}
H(\omega) &= \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \\
&= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n}
\end{align}
$$

To simplify this, we can use the known formula for the sum of the first $N$ terms of a geometric series:

$$
\sum_{n=0}^{N-1} e^{-j\omega n} = \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}}
$$

What we're most interested in for your case is the magnitude response of the filter, $|H(\omega)|$. Using a couple simple manipulations, we can get that in an easier-to-comprehend form:

$$
\begin{align}
H(\omega) &= \frac{1}{N}\sum_{n=0}^{N-1} e^{-j\omega n} \\
&= \frac{1}{N} \frac{1-e^{-j \omega N}}{1 - e^{-j\omega}} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{e^{j\omega N/2} - e^{-j\omega N/2}}{e^{j\omega /2} - e^{-j\omega /2}}
\end{align}
$$

This may not look any easier to understand. However, due to Euler's identity, recall that:

$$
\sin(\omega) = \frac{e^{j\omega} - e^{-j\omega}}{j2}
$$

Therefore, we can write the above as:

$$
\begin{align}
H(\omega) &= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{j2 \sin\left(\frac{\omega N}{2}\right)}{j2 \sin\left(\frac{\omega}{2}\right)} \\
&= \frac{1}{N} \frac{e^{-j \omega N/2}}{e^{-j \omega/2}} \frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}
\end{align}
$$

As I stated before, what you're really concerned about is the magnitude of the frequency response. So, we can take the magnitude of the above to simplify it further:

$$
|H(\omega)| = \frac{1}{N} \left|\frac{\sin\left(\frac{\omega N}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}\right|
$$

Note: We are able to drop the exponential terms out because they don't affect the magnitude of the result; $|e^{j\omega}| = 1$ for all values of $\omega$. Since $|xy| = |x||y|$ for any two finite complex numbers $x$ and $y$, we can conclude that the presence of the exponential terms don't affect the overall magnitude response (instead, they affect the system's phase response).

The resulting function inside the magnitude brackets is a form of a Dirichlet kernel. It is sometimes called a periodic sinc function, because it resembles the sinc function somewhat in appearance, but is periodic instead.

Anyway, since the definition of cutoff frequency is somewhat underspecified (-3 dB point? -6 dB point? first sidelobe null?), you can use the above equation to solve for whatever you need. Specifically, you can do the following:

Set $|H(\omega)|$ to the value corresponding to the filter response that you want at the cutoff frequency.
Set $\omega$ equal to the cutoff frequency. To map a continuous-time frequency to the discrete-time domain, remember that $\omega = 2\pi \frac{f}{f_s}$, where $f_s$ is your sample rate.
Find the value of $N$ that gives you the best agreement between the left and right hand sides of the equation. That should be the length of your moving average.


Standard Deviation in Gaussian Blur I have a function that performs gaussian blur on image for some specific $\sigma$ (the standard deviation).

It first computes kernel of size $\lceil 3\sigma \rceil$ and then performs convolution with that kernel.

However, I would like to specify blurring radius in pixels rather than $\sigma$.

I suppose that the blur radius (in pixels) is just $\sigma^{2}$ as this denotes variance of a random variable.

Is that right? Can the same thought be extended to 2D?

UPDATE:

The problem is that I need to do things like building a gaussian pyramid (successively blurred and downsampled image).

When the image gets downsampled to 1/2 of its width, I suppose I need a gaussian blur of radius 2 pixels ($\sigma=\sqrt{2}$ ?). And for 1/4 subsampling, I would need blur of 4 pixels ($\sigma=2 ?)$... But I am not sure about that...

Yes, you can always use pixel units for $\sigma$. In 2-D you will have $\sigma_X$ and $\sigma_Y$.


Frequency Spectrum with c++ I have written a script to plot the frequency spectrum of a sinusoidal signal with c++. Here are the steps
1-Applying Hanning window
2-Apply FFT using fftw3 library
I have three graphs: Signal, Signal when is multiplied to Hanning function, and the frequency spectrum. The frequency spectrum looks wrong. It should have a pick at 50 Hz.Any suggestion would be appreciated. Here is the code:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

int main()
{
int i;
double y;
int N=50;
double Fs=1000;//sampling frequency
double  T=1/Fs;//sample time
double f=50;//frequency
double *in;
fftw_complex *out;
double t[N];//time vector
double ff[N];
fftw_plan plan_forward;

in = (double*) fftw_malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

 for (int i=0; i< N;i++)
 {
    t[i]=i*T;
    ff[i]=1/t[i];
    in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
    double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
    in[i] = multiplier * in[i];
  }

  plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );

  fftw_execute ( plan_forward );

  double v[N];

  for (int i = 0; i < N; i++)
    {

    v[i]=20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])/N/2);//Here I have calculated the y axis of the spectrum in dB

    }

   fstream myfile;

   myfile.open("example2.txt",fstream::out);

   myfile << "plot '-' using 1:2" << std::endl;

   for(i = 0; i < N; ++i)

    {

      myfile << ff[i]<< " " << v[i]<< std::endl;

    }

 myfile.close();

 fftw_destroy_plan ( plan_forward );
 fftw_free ( in );
 fftw_free ( out );
 return 0;
  }

I have to add that I have plotted the graphs using gnuplot after inserting the results into example2.txt. So ff[i] vs v[i] should give me the frequency spectrum.

Here are the plots:   Frequency Spectrum and Sinusoidal time Window respectively:

Edit: After applying the changes to frequency sampling, this is the graph I got:
Edit: The first graph is the plot for frequency spectrum vs time

You seem to be doing something odd with your frequency axis. This:

t[i]=i*T;
ff[i]=1/t[i];

should probably read:

t[i]=i*T;
ff[i]=Fs*i/N;

The frequency axis you have now is nonlinear. The FFT generally generates linearly spaced frequency samples.


How can I find the transform function, amplitude and frequency response of a digital filter in python? I have applied a laplace filter mask to an image and now I want to find the amplitude and freqency response of a laplacian filter: [[1,1,1], [1,-8,1], [1,1,1] ]. I know I need to first find the transfer function, however, I am unable to do this as well programatically in python.

Here is the code I have currently:

from matplotlib.pyplot import imread
import matplotlib.pyplot as plt
from scipy.ndimage.filters import convolve
import numpy as np
import cv2 as cv

laplaceMask = np.array([ [1,1,1], [1,-8,1], [1,1,1] ])
img = imread("download.jpg")

channels = []
for channel in range(3):
    result = convolve(img[:,:,channel], laplaceMask)
    channels.append(result)

img = np.dstack((channels[0], channels[1], channels[2]))

plt.imshow(img)
plt.show()

Assuming each element is the filter coefficients for the R, G, and B channels (I do not know if this is indeed the case here), you can get the frequency response of each channel using scipy.signal.freqz:

scipy.signal.freqz(c, 1, v)

Where c is the filter coefficients and v is the signal data for that channel.


Can you quarter the processing time for real, symmetric FFTs? The DFT of a real signal is Hermite-symmetric, so you can roughly halve the computation time/memory by not bothering to calculate half the values of the spectrum (and complex conjugating the existing values and copying them to the second half if needed).  So the rfft operation takes N samples and outputs N/2 spectrum bins in half the time, for instance.

Signals which are even-symmetrical and real have even-symmetrical and real spectra (and real spectra take half the memory to store as complex spectra), so for symmetrical input, can the calculation be halved again by only using the first half of the signal (N/2) to generate the first half of the spectrum (N/2)?  How?

Is there a way to take the regular FFT of the N/2 half-signal and manipulate the output to produce the N/2 half-spectrum?

(real/odd ⇔ imaginary/odd would work, too, but real/even ⇔ real/even case is simpler to follow.)

I'm not sure about quartering the process time but this is what I can think of for now. N point DFT of x.

$X_k = \sum_{n=0}^{N-1}x_nW_N^{kn}$

By the even symmetry assumption

$X_k = \sum_{n=0}^{N/2-1}x_n(W_N^{kn}+W_N^{k(N-1-n)})$

Then we simplify

$X_k = W_N^{k(N-1)/2}\sum_{n=0}^{N/2-1}x_n(W_N^{kn}W_N^{-k(N-1)/2}+W_N^{-kn}W_N^{k(N-1)/2})$

$X_k = 2W_N^{k(N-1)/2}\sum_{n=0}^{N/2-1}x_n\cos(2\pi k\frac{n-(N-1)/2}{N})$

There are probably some errors but if simplifications are valid the 'FFT' size is reduced by a factor of 2 and the twiddle factors are real (this maybe corresponds to your quartering). As in normal FFT the symmetry and periodicity can be exploited to decompose the 'symmetric FFT' into radix-2 steps for instnace.


Vibration Analysis by using FFT I have 6 column vectors $(\mathbf A_x,\mathbf A_y, \mathbf A_z, \mathbf G_x, \mathbf G_y, \mathbf G_z)$ taken by accelerometer to analyze the vibration of a construction. Each column is a $64998\times 1$ valued-vector and overall time is $664$ seconds. I've done something but I don't know how I can fix it.

Fs=97.885;
L=64998;
NFFT=2^nextpow2(64998);

AXdft = fft(Ax);
AXdft=AXdft(1:NFFT/2+1,:);

freq = linspace(0,Fs,NFFT/2+1);

plot(freq,2*abs(AXdft));
xlabel('Hz');ylabel('Magnitude');

figure;
plot(freq,20*log10(abs(AXdft)));
xlabel('Hz');ylabel('dB');
title('Plot in dB');

BTW I will add each column over other via hold on

getting a figure but missing its size. The frequency doesn't appear in figure and I guess I must do PSD, Can anyone help?

I believe that your problem is using nextpow2 to set the fft length.  fft(x) returns a vector length(x). fft(x) is not constrained to a power of 2. If you want to set a particular size by zero padding, use fft(x,N). NFFT=length(Ax) is the actual size of DFT.

Using psd has advantages but if you use the proper size, your code will work (assuming I am correct) . Taking a single long DFT isn't an accepted way to estimate a psd, unless you are looking for very fine, high snr features.


Does anyone know how and when and where used Legendre polynomial? I'd like approach to Legendre polynomial.
But i don't make concept with Legendre polynomial.
So if you know where and when does L.P used and how L.P used, would you please let me know about this? To more understand it.

Hmm.. specifically in the application.

Legendre polynomials are orthogonal on the interval $[-1,1]$ and are therefore useful for approximating functions on a finite interval in a least squares sense. I.e., if you are looking for a polynomial $p^*_K(t)$ of degree $K$ which approximates a given function $f(t)$ on the interval $[-1,1]$ - other intervals are possible by shifting and scaling - by minimizing the error

$$\epsilon=\int_{-1}^1[f(t)-p_K(t)]^2dt$$

then the optimal polynomial $p^*_K(t)$ is given by

$$p^*_K(t)=\sum_{k=0}^Ka_k\bar{L}_k(t)\tag{1}$$

where

$$\bar{L}_k(t)=\sqrt{\frac{2k+1}{2}}L_k(t),\quad k\in\mathbf{N}$$

are the normalized Legendre polynomials $L_k(t)$ (i.e., $\bar{L}_k(t)$ are orthonormal, not only orthogonal). Because of the orthonormality of $\bar{L}_k(t)$ the coefficients $a_k$ in (1) can be easily found by computing

$$a_k=\int_{-1}^1\bar{L}_k(t)x(t)dt$$


SNR calculation when noiseless signal is not known I have a noisy signal. Then using filtering techniques I got the estimated noiseless signal. I dont have the original noiseless signal. I have calculated SNR the following way:
$$
\textrm{SNR}=10\log\left(\frac {\textrm{Estimated Noiseless Signal}^2} {\left(\textrm{Estimated Noiseless Signal-Noisy Signal}\right)^2}\right)
$$

Is this method right?

Assuming you're calculating the $\textrm{SNR}$ in decibels, it is defined by the following equation:

$$
\textrm{SNR}_{db} = 10\log_{10}\left(\frac{P_{signal}}{P_{noise}}\right)^2
$$

You currently have (with squares consolidated), where $P$ is the noisy signal:

$$
\textrm{SNR}_{db} = 10\log_{10}\left(\frac{P_{signal}}{P_{signal}-P}\right)^2
$$

This doesn't look right because you will end up with a negative denominator since $P>P_{signal}$, where in essence $P=P_{signal}+P_{noise}$. Although it ends up being the same thing because the squaring forces it to be positive, it seems a bit of a fudge.

I think what you meant to do was to subtract the estimated signal from the noisy signal to estimate the noise power:

$$
\textrm{SNR}_{db} = 10\log_{10}\left(\frac{P_{signal}}{P-P_{signal}}\right)^2
$$

One way of obtaining your $P_{signal}$ is to understand what your transmitter is capable of putting into the medium. It would then follow that the further away your receiver was to your transmitter, the lower the $\textrm{SNR}$. Remember, a low SNR can mean two things, either your transmitter isn't loud enough, or the background noise is too great.


SINR averaging in communication system Under the simulation of a communication system suffering from noise (AWGN) and co-channel interference, what is the appropriate way of calculating SINR?
Could one calculate it at a packet level as:

$$
SINR_{packet} = \frac{\sum{\left \| x \right \|^{2}}}{\sum \left \| n \right \|^{2}+\sum \left \| z \right \|^{2}}
$$

where $x$ is the total useful waveform for that packet, $n$ is the corresponding noise waveform and $z$ is the waveform of the interferer, or is it necessary to calculate it on a per bit basis as

$$
SINR_{bit} = \frac{\sum{\left \| x_{b} \right \|^{2}}}{\sum \left \| n_{b} \right \|^{2}+\sum \left \| z_{b} \right \|^{2}}
$$

where $x_{b}$ is the useful waveform for that particular bit, $n_{b}$ is the corresponding noise waveform for that bit and $z_{b}$ is the waveform of the interferer for that bit and then take the average across all bits of the same packet like so:

$$
SINR_{packet} = \frac{1}{k}\sum_{bit=1}^{bit=k}{SINR_{bit} = \frac{1}{k}\sum_{bit=1}^{bit=k}{\frac{\sum{\left \| x_{bit} \right \|^{2}}}{\sum \left \| n_{bit} \right \|^{2}+\sum \left \| z_{bit} \right \|^{2}}}}
$$

Which one is supposed to converge more quickly and give more realistic results? Which one would be better related to bit error rate (BER) and which to packet error rate (PER)? Finally could Central Limit Theorem be of any use proving that the first and last expression converge to the same number for large number of bits in a packet?

I believe real systems calculate SNR as in the first formula, except of course for the $z$ term.

Signal-to-interference-plus-noise ratio (SINR) is by definition

$SINR = \frac{P}{I+N}$,

where $P$ - is useful signal power, $I$ - interference signal power and $N$ - noise power in the band of interest. Look

http://en.wikipedia.org/wiki/Signal-to-interference-plus-noise_ratio

for example.

So the first formula you've presented seems to be rigth one. Time span of calculation $SINR$ depends on your needs and interference nature. The common way is using duration of $x$. In general you have to distinguish both $x$ and $i$ at first, because $i$ is non-stationary random process so you can't think it as receiver noise. So in most common case in real life you can estimate $SINR$ only if you are able to combat such an interference in you system.

If inteference appears as short pulses in time you can do measurement with relatively short sliding window to estimate $SINR$ in dynamics. If signal power is quite stationary in short averaging you can use previously estimated $P$ in formula above to find out $SINR$ even if demodulator fails to defeat it.

You can combine measurement of $SINR$ at the input of demodulator, $MSE$ at its output and $BER$ at the decoder's output to make some sophisticated modem quality control.


Simulating noise in computed tomography reconstruction In the research of the computed tomography (CT) reconstruction, one needs to simulate the noise during CT projection and capturing. Then, my questions are:

How about the proper noise type?
How to add the noise, in the sinogram domain or the frequency domain?

Any reference materials are also preferred.

The easiest, and probably most straight forward way, seems to add the noise in the measurement domain, hence the sinogram.


OpenCV detectin skin colour invariant with illumination changes I have to capture skin colour independent of illumination.
I used the following criteria(using hsv and rgb) but it ain't working:

int h = get_hue(...);
int s = get_saturation(...);
int r = get_red_component(...);
int g = get_green_component(...);
int b = get_blue_component(...);

if ((h<38)&&(h>6)&&(s>0.23)&&(s<0.68)){
    // skin color
}
else if ((r>95) && (g>40) && (b>20) &&
         ((maximum(r,g,b)-minimum(r,g,b))>15) &&
         (abs(r-g)>15)&&(r>g)&&(r>b)){
   // also skin color
}

It works in the morning when my skin is illuminated but by evening it does not work.

Any Help will be appreciated. Thanks in Advance.
(PS -my skin ain't white.)

In my experience best method for this is converting it to Lab color space.
L represents the light, and a and b are light independent.
OpenCV supports Lab color scale conversion.


What is filter memory? I have just started learning some DSP, I keep seeing things like memory in regard to filters.  I am unsure exactly what this means and why it is useful.  I understand that increasing the number of taps increases the memory etc.. and a tap is essentially just a coefficient for changing the signal.

Filter memory refers to the effect of past input on the current input. Infinite impulse response (IIR) filters and finite impulse response (FIR) filters are fundamentally different from the perspective of their memory. The memory of an FIR will essentially be the number of taps, while IIR has theoretically infinite memory (although when they are stable, the memory passes under the numerical accuracy of the system beyond a certain point).

More on this http://gregstanleyandassociates.com/whitepapers/FaultDiagnosis/Filtering/Filter-Memory/filter-memory.htm.


Radar type application for narrow band only? Not sure if this is the proper forum to ask this question, but still I am trying my luck.
I am looking for an application (in any area, medical, defence, or any other area) if exist, such that we want to dectect something by sending a waveform and get a retun on a linear array and determine the exact location of the target from the return signal. The catch is due to the situation or the physical limitation we cannot use wideband waveform. (We know that linear FM or wavefrom with higher bandwidth will help in spatial resolution, but we cannot use it due to some constraint.) In addition, the sampling rate cannot be increased so that we cannot form correlation matrix for better spatial resolution. This may also mean that target changes fast and we must rely only on time snapshot.
It appears in such case we stuck with Fourier beamforming.
The question is: Is there any such application where the above constaint really exists.
Any ideas will be highly appreciated.

These constraints absolutely exist. There are the norm! We could only wish in our wildest dreams to use as wide as bandwidth as we like. There are many areas in a radar system that place limitations on how wide the bandwidth can be and we'll go over a few straight forward ones. Mainly we're talking about limitations due to the antenna and waveguide as well as signal processing requirements.
Antenna and Waveguide
With antennas, especially arrays, their geometries can be defined in terms of wavelengths. Given that wavelength, and therefore frequency, we can yield the antenna's directivity pattern. It's important to stress that when you typically look at an antenna pattern, it is for a single RF frequency only.
Even simple pulses have some bandwidth around the RF frequency but are usually not too large. This allows us to make the narrowband approximation that the antenna pattern is the same for frequencies near the RF frequency, which is how the pattern was defined.
With an exceptionally wideband waveform the frequencies transmitted are now farther apart and so the antenna pattern itself can change significantly intra-pulse. These affects are more extreme at lower RF frequencies where desired bandwidths are well within the order of the transmitted frequency.
These effects are undesirable because (including but not limited to):

The effective radiated power is no longer the same during the pulse.
Because the antenna pattern itself is changing during the pulse, it distorts the waveform itself and you incur SNR degradation as well as undesired phase behavior which may lead to other issues, such as errors in angle measurements.

The waveguide itself, as well as other receiver/exciter components, need to support the bandwidth of interest and thus will require the use of certain technologies and processes to support it. This is a nice way of saying that you're going to be spending a lot of money and making sacrifices in size, weight, and power (SWaP).
Signal Processing
Assuming that your RF hardware is taken care of, there are issues with data collection and front-end signal processing. High bandwidth waveforms require high sampling rates. Introducing higher sampling rates has the following overarching issues

High sampling rate analog-to-digital converters (ADCs) can be prohibitively expensive...for a good one.
Depending on the bandwidths and pulse widths desired, you may have to deal with a large number of samples. In critical applications such as a fire control radar your computing resource (usually FPGAs/CPUs) need to be able to take care of common tasks such as pulse compression, FFTs, memory copies, etc. within a very strict time schedule.

These are just a few examples. I hope it's clear that the narrowband restriction is many times the norm when designing a system. Thankfully, technology is advancing where we can get away with higher bandwidths easier than we could in the past.


If we define samples manually how could there be a Sampling Frequency? I just got a thought when generating matrices , plotting it and taking its FFT,

Lets say, I have randomly input a matrix of length 256,

[1 7 9 5 7 1 -2 0..................8]

According to the FFT, the resulting 256 FFT coefficients correspond to the following frequencies:

0, $f_s/256$, $2 f_s/256$, $3 f_s/256$, $\ldots$, $255 f_s/256$

but since the samples i taken has no sampling frequency defined how could there be some sampling frequency ? Like in an image ,

Since the image is just a 2 dimensional Matrix , if i randomly generate an image in a matlab with random pixel values (sample points) and plot its FFT , the matlab will show the frequencies of that 2-D matrix but since there is no Sampling Frequency is defined how matlab calculates the frequency for each of the coefficient ?

Units can be arbitrary, such as ISO Meters, or the length of the current king's foot (maybe longer than last year after the assassination of some shorter king).

A good default unit if there isn't some ISO standard, useful natural unit, or you can't think of one to make up (how long is your monitor's display, or your left foot, etc.), is the length of your window of data.  e.g. 1 cycle per your vector of 256, 2 cycles, 3 cycles, etc. per your vector length.  If your vector displays as 3 inches long on your display, then you can also easily convert cycles/vec into cycles/inch, or per/mm, per/fortnight, etc.


Why I get another spectrum when upconvert an o-qpsk baseband signal with a bladerf compared to a hackrf device? I have just start working with bladerf instead of former used hackrf hardware and I need to solve a problem. Currently I use the following bladerf hardware, libraries and firmware:
bladerf (2.0 micro xA9)
libbladeRF version:         2.2.0-2018.12-rc3-2-ppabionic
Firmware version:           2.3.2
FPGA version:               0.10.2 (configured by USB host)
I tried to simply upconvert an oqpsk baseband at 1100 MHz (4 Msps) with bladerf device using
gnuradio (GRC 3.711) block osmocom-sink and I got a unexpected and pretty wierd baseband spectrum compared to hackrf hardware.
For better understanding, I have attached some hardcopies:

From my point of view, I made identical setups for hackrf and bladerf device on osmocom-sink block. Later on downconverting and demodulation works very well with hackrf hardware, with bladerf device it totaly fails.
Any suggestions, I need to check or I just made wrong ? May I need to calibrate the bladerf hardware first for proper baseband processing ?
Many thanks in advance.
Best regards
calenberger

One source of this is spectral regrowth (overdriving a device into saturation such that the pulse shaping that limits the spectrum is degraded), but compared to typical spectral regrowth spectrums I have seen, this looks more like your signal was greatly attenuated and then amplified together with noise, and then that combined signal passed through a bandpass filter. Check if there is a 5 MHz filter earlier in the Tx chain to further support this theory, and then also try to reduce the power level at the input to the device.


What is the purpose of a window function? I've read that a window function is a mathematical function that is zero-valued outside of some chosen interval.

Can someone please explain what this really means?

I've read that a window function is a mathematical function that is zero-valued outside of some chosen interval.
Can someone please explain what this really means?

It means exactly what it says, but I assume you really want to know why do we care?
The purpose of window functions in signal processing is to remove (hopefully) extraneous information, so that subsequent analysis produces more meaningful results. Checkout wikipedia (and other sources) for specifics, but windows functions do pre-processioning in the time domain to help improve results of later processing in the frequency domain.
This hopefully helps explain why there are so many different window functions.  Basically, different window functions are designed to help with different problems present in the time domain data.  Often there is no right/wrong window functions, there are just more/less useful functions for solving the problem at hand.


Internet from Direct Broadcast Serivce Internet is very much costly in India , like 3-5 USD  per GB @ 3G networks and even in broadband 100 USD for a unlimited banwith for 2 MBps but here DTH (Direct-To-Home broadcast service ) is very common and cheap , like 4$ for a month with most of the channels . I think to buffer a 480p video we would need 200 KB/s (2MBps ) plan in internet but how DTH providers give that much bandwidth for TV service for such a low price .

My main point here is :

Can we use DTH service for Internet using (direct from satellite like the DTH works).
Can we get same bandwidth for using internet just at same price just in case of DTH .

No you can't, because the internet depends on you having a method of sending data. For example, the question you've asked here came here somehow, so it's obvious you need to be able to send data.

Now, TV is broadcasting, it's a one-way transmission, without a back channel.

There are satellite internet access services, but they are far more costly than what you've described, and much much slower.

Also, some countries had things that used a telephone line as uplink and a sattelite downlink, but they suffered from cost, the asymmetry between the fast satellite downlink and the slow telephone uplink, and the poor latency.


Multiplication in the wavelet domain, what does it look like in real space? Say you wavelet transform a signal, multiply by some values, and untransform. Is there anything we can say about the effect in real space? I know that it's not exactly a convolution, because the convolution theorem doesn't really exist for the wavelet transform. So if it's not a convolution, what is it?

Also, is the effect different if you multiply in an orthogonal wavelet basis (at, say, dyadic frequencies) versus a continuous transform at all frequencies?

They can be  used for time-varying filtering of non-stationary signals, under names like Gabor multipliers, STFT filters for instance.  You can find a review reference in the chapter Linear Time-Frequency Filters, 11.1.3 Implicit design and generalization in Spreading function representation of operators and Gabor multiplier approximation.

You also have the 2004 paper The convolution theorem for the continuous wavelet transform. My  experience in this domain is limited to localized matched filtering, by  multiplying coefficients by scalar complex numbers to locally modify amplitude and phase (with complex continuous wavelet frames), see Adaptive multiple subtraction with wavelet-based complex unary Wiener filters.

The redundancy of the continuous transform certainly helps a lot. With DWT (a variant of them) I have also used  masks to remove localized high amplitude disturbances in some discrete wavelet subbands. The masks were initially binary ($0$ on the places to remove, $1$ outside), but they were tapered to have a smooth profile between $0$ and $1$. Not very theoretical, but efficient to separate directional data.

Finally, you can find across-scale multiplication for ROI (region-of interest) coding with JPEG 2000.


Is it only me or matlab's periodogram's function is confusing? I am really confused with matlab's periodogram.

First of all, it says in the documentation that its unit is "dB per unit frequency". What does that really mean? I just don't quite get it.

Second of all, I'm trying to evaluate an ADC. I tried to measure its linearity but that won't work because when I reduce the level of the input signal by ,let's say, -3dB I excpect to see -3dB reduction in the highest peak of the spectrum, but that's not the case.

Last but not least, with this "dB per unit frequency" when is the highest peak actually 0 dB? It should reach that at a point, right? Because as far as I know dB is always a ratio of some kind. How can I modify it as to give me the values in dBFs for example? so I could get a sense of what my data look like.

Sorry for the long post, Mathworks' Periodogram just doesn't say it clear.

Thank you all.

As the documentation states, periodogram provides a power spectral density estimate pxx:

[pxx, w] = periodogram(x);

meaning that it shows how the total variance of the signal, var(x), is distributed over the frequency w. This implies that the integral of the estimated spectral density over the frequency axis

trapz(w, pxx)

should correspond to that total variance.

This also means that the psd is not normalized, neither to a total variance of 1, nor to a maximum peak of 1 (= 0 dB). That the result is plotted (not computed!) in dB merely means that the vertical axis shows $10 \, \log_{10}($ pxx $)$, not that it shows a ratio. There is therefore no reason to expect a peak to align with 0 dB; depending on the range of values of the signal it can be both much higher and much lower.

To produce a plot from pxx and w, use

plot(w, 10 * log10(pxx))
xlabel('normalized frequency   [rad / sample]')
ylabel('power spectral density   [(rad/sample)^{-1}] in dB')

Reducing the level of the total signal by 3 dB means to divide the signal x by a factor of 2. If you do this numerically

periodogram(x)
hold all
periodogram(x / 2)

you will see that the second psd estimate has exactly the same shape as the first, it's just 3 dB lower. If this is not what you find with your measurement data, then this means that either your volume regulation or your ADC do not work correctly.

If you are really sure you want a normalized spectrum, do the following: First, compute the psd estimate:

[pxx, w] = periodogram(x);

Then, normalize the psd, e.g. with respect to the highest peak:

npxx = pxx / max(pxx);

You can then plot the normalized psd on a logarithmic scale like this:

plot(w, 10 * log10(npxx))
xlabel('normalized frequency   [rad / sample]')
ylabel('normalized power spectral density')

Note that npxx does not have a unit.

A note about units: The unit of numbers returned as pxx by periodogram & friends depends on the units of the input signal  x. Let's say these are dimensionless numbers; then the unit of the numbers in the pxx output is $\rm \frac{1}{rad / sample}$, because it is a density over w which is in units of $\rm rad / sample$. Or let's say these are voltage measurements coming from a calibrated ADC, so that the numbers in x are in volts. Then the unit of the numbers in the pxx output is $\rm \frac{V^2}{rad / sample}$.

If you additionally specify the sampling frequency in the call to periodogram:

[pxx, f] = periodogram(x, [], [], fs);

and fs is in Hz, the pxx is a density over frequency f, and its unit is $\rm \frac{V^2}{Hz}$.

Now if pxx is plotted, it is often plotted logarithmically, like stated above. If the function $10 \log_{10}$ is used, then it is said that the result is in decibel, dB. Then sometimes the unit of 10 * log10(pxx) is written as $\rm \frac{dB}{rad / sample}$. Strictly speaking, this expression doesn't make sense because dB is not just a unit (a different measurement stick) but denotes the use of a different quantity: the logarithm. That's why above I wrote the unit as [(rad/sample)^{-1}] in dB.


What are the advantages and disadvantages of Kalman filter compared with FIR, IIR and low pass filter to filter data with noise? It is known that the Kalman filter can filter the data with noise. I also find it works well after using it compared with FIR, low pass filter,etc. Now, I have a couple of questions about the advantages and disadvantages of the Kalman filter compared with FIR, low pass filter,etc. (In accordance with the conditions of using Kalman filter)

Can we deem that traditional filters such as FIR and low-pass filter are designed to be used at specific frequencies, while kalman filter can filter signals with noises at any frequency as long as the computing speed is large enough and the system model does not change?
The kalman filter can adjust the Kalman gain according to the actual measurement accuracy, so as to obtain the optimal solution. Is it right to think that Kalman filter is more intelligent than FIR?
When the measuring instrument itself has errors, is it right to think the data obtained by using kalman filter is more accurate for the Kalman filter can adjust the weights of system model and measure respectively?
When the noise of measurement is large (the instrument is not accurate), can the adoption of Kalman filter (the model is well established) obtain more accurate data than FIR, etc. That is, can the kalman filter be used to make up for the error caused by the low precision of the measuring instrument by establishing an accurate model?
The core of kalman filter is a set of iterative equations. Can we think the kalman filter has the advantages of simple design and strong universality compared with FIR? Of course, it is based on the establishment of a good system model and measurement model.

What I understand is not necessarily right or sufficient. Please help me to answer my questions or give me any suggestions. Thanks!

Kalman filters really aren't that special, and you seem to be missing the point of  a Kalman filter.  A Kalman filter is really just a generally time-varying, generally IIR, generally multi-input multi-output filter that's been designed using a specific procedure.

Can we deem that traditional filters such as FIR and low-pass filter are designed to be used at specific frequencies, while kalman filter can filter signals with noises at any frequency as long as the computing speed is large enough and the system model does not change?

Yes and no.  Yes, typically time-invariant IIR and FIR filters are designed starting with a specification in the frequency domain -- but they don't have to be.
A Kalman filter is simply a filter that's designed based on a model of a process, to be optimal under a certain set of rules.  But you can give me a filter specification for an IIR or FIR filter, and I can define a system model that will cause you to design a "Kalman filter" that asymptotically settles out to being exactly an IIR or FIR filter with that specification.  Then if you make your "Kalman filter" into a steady-state Kalman, it will be that original IIR or FIR filter -- just with more (and more obscure) work behind it.

The kalman filter can adjust the Kalman gain according to the actual measurement accuracy, so as to obtain the optimal solution.

Well, no.  The Kalman filter can adjust the gains according to the estimated or understood measurement accuracy.  Depending on your religion, either no one, or God only, knows the actual measurement accuracy.

Is it right to think that Kalman filter is more intelligent than FIR?

I have never had an intelligent conversation with either.  Zero equals zero, so no.  A Kalman filter may be more applicable than a FIR filter in certain circumstances -- but often it is not.

When the measuring instrument itself has errors, is it right to think the data obtained by using kalman filter is more accurate for the Kalman filter can adjust the weights of system model and measure respectively?

Not necessarily.  Particularly since Kalman filters have a well-known tendency to lack robustness in the face of inaccuracies of the system model -- get the system model wrong (and you will, because you're not God), and the filter won't be robust.
Sometimes it's far better to use a simple filter that's good enough, than a super-duper-fancy filter that would be really good if only you'd designed it right, based on information you don't have.

When the noise of measurement is large (the instrument is not accurate), can the adoption of Kalman filter (the model is well established) obtain more accurate data than FIR, etc. That is, can the kalman filter be used to make up for the error caused by the low precision of the measuring instrument by establishing an accurate model?

"Kalman" is a Hungarian surname that people whose ancestors come from a certain area of Germany tend to have.  It does not mean "magic".  A Kalman filter is just one way of making an optimal filter -- but it cannot do better than optimal.
Particularly for single-input, single-output systems, a plain old IIR or FIR filter, specified in the frequency domain, can do a better job (because of robustness issues) than some laboriously designed Kalman.

The core of kalman filter is a set of iterative equations. Can we think the kalman filter has the advantages of simple design and strong universality compared with FIR? Of course, it is based on the establishment of a good system model and measurement model.

We can think of the Kalman filter as being a right pain in the behind to design correctly, and one that forever more demands that you have staff on hand that understands it before your product can be adequately maintained.  If simple will work, simple should be adopted -- and based on the number of people who can understand IIR & FIR filters vs. the number of people who truly understand how a Kalman filter works, I'd say that a Kalman filter is not "simple".


STFT - Putting it all together I'm currently reading this paper: http://www.rssd.esa.int/SP/LISAPATHFINDER/docs/Data_Analysis/GH_FFT.pdf
On page 21, it has a section called "Putting it all together", it states the following steps:

Input data
Frequency resolution and length of DFT
Window function
Splitting of the data stream
FFT
Scaling

Now, the "Window function" is mapped out before the "Splitting of the data stream" if this is the case, should I therefore multiply each of the bins within the data stream by each of the window values (Hanning)?
For example, at the moment, my Hanning window size is: 256 and the data is split into segments of size 256x128 in this case.. N = size of stream and then compute hanning(N)
I've seem a lot of DSP engineers split the window up into frames, before applying the window?
Thanks

If your window is a rectangular window, then there is no mapping out since splitting the stream into frames is effectivley applying that window. Once you split the stream into frames, if your window is not rectangular, you would then multiply the window against the frame.

y[i] = x[i] * w[i], for i = 0 to N-1, then take the FFT of y.


Simulating musical instruments using their impulse response I studying DSP and trying to find the limitation of 'teaching' a software to emulate a musical instrument by sampling its notes at different velocities and building an LTI system for simulating it.

Is that feasible? is that even LTI? is that limited to certain instruments only? what are the problems I'm about to encounter trying to do that?

Thank you

there are a lotta people doing research regarding music synthesis and DSP.  yes you can sample notes and different key velocities and play them back with different MIDI velocity parameters.  if you wanted to somehow interpolate between a note played at $mf$ to $f$ and on to $ff$ and $fff$urther, you would need to find a way to phase align the waveforms.  IMO the only decent analytical way to do that in the time domain is wavetable synthesis where the wavetables can be lined up (or it's frequency domain, it's called additive synthesis) - additive synthesis is more general, wavetable requires the overtones of the note to be mostly harmonic.  doing this stuff right is not particularly easy.

on the other hand, since you brought up "impulse response", there are folks who do physical modeling of musical instruments, where they model the instrument as a form of resonant LTI filter and they separate the impulse response of the LTI filter from the driving function, which is not an impulse but something more.  then, for different key velocity hits, they have different driving functions but the same LTI impulse response.  sometimes it's called commutated synthesis or https://ccrma.stanford.edu/~jos/pdf/svd95.pdf.


How to calculate the displacement from the measured acceleration data? I have an acceleration data which is collected from a vehicle. I do also have a displacement measurement.

I would like to calcuate the displacement vector from the measured acceleration data, but no success yet.

Here is the data (the previous link was not working I updated the link):
https://ufile.io/zyw3r

What I try:

filename = ('acc_time.csv');
accdata = csvread('acc_time.csv');

[mydata, myheader] = xlsread(filename);
for i = 1:length(myheader)
      % compose a command to assign each column to a variable with the same
      % name as the header
      commandExec = [myheader{i}, ' = ', 'mydata(:,', num2str(i) , ');'];
      % execute the composed command to actually create the variable
      evalin('base', commandExec );
end

vel_vector = cumtrapz(time,acc);
disp_vector= cumtrapz(time,vel_vector);
plot(time, disp_vector)

and I get:

but the displacement vector should be squiggy not like linear incremental.

Where am I wrong?

Generally, the observed effect occurs when there is an offset and/or linear drift present in your measured data. Double integration then leads to the quadratic (or higher order) effect observed in the data.

As @Stanley Pawlukiewicz states in his answer this can be due to the accelerometer registering gravity (not all accelerometers are capable of registering a constant acceleration field) but also due to temperature fluctuations, (electro)magnetic interference, strain at the base of the accelerometer, etc. Sometimes, it’s possible to “shield” your sensor for these effects but this can also be extremely challenging or expensive. In the latter case pre-processing your data can be an option too enable extracting meaningful results from your data.

The most straightforward approach to get rid of the quadratic drift in your double integrated acceleration signal is to detrend the data before integrating. If you substract the mean value before integrating the acceleration and velocity signals you should get reasonably “close” to the desired result (see matlab’s detrend function for other options).

Note that there also appears to be other noise and/or content in your measured acceleration signal; i.e. you will probably be able to get “close” to the desired result but will not get the exact same position signal. If that is critical you could use a more sophisticated signal processing approach or, if possible (and preferably) a better sensor (e.g. a position sensor).


Why choose anything but highest compression for PNG? I just installed LibreOffice for the first time. In the Impress
application for presentation slides, when I choose to save an image on
a slide, the file type is PNG and there is a slider for compression,
positioned approximately midway at 6. Internet searching reveals that
PNG compression is nonlossy.  Why would anyone choose anything but the
highest compression?

Simply because the highest compression typically is significantly more CPU-intense (it tries out multiple different approaches to represent successive lines).

This really shouldn't make much difference on a modern PC for saving a few images. Then again, in practice, libpng seems to be pretty slow, so this might make a difference, especially for people working on older devices.


How to apply dictionary learning to several images? I'm implementing a dictionary learning algorithm for doing sparse coding on images and I have a problem how to learn a dictionary corresponding to several iamges.

Let’s say X is the data matrix. If I apply dictionary learning only to one image then X is equal to the image. But what if I want to learn the dictionary for several images ? Should I in this case just put all images into X, i.e. X = [image1, image2, image3…]?

Another approach I thought of was to learn a dictionary for each image separately and then decide online which dictionary to use (by e.g. comparing histograms or some similarity measures to decide to which training image the new input image is most similar to and then choosing the corresponding dictionary).

What appraoch should I use?

Dictionary learning is usually done on more than one image. So yes, I would normalize and concatenate multiple images together. What I usually do is create 8x8 patches from images and vectorize them so that they represent one column (data-item) from the training set (X in your case).


Beginner LTI system Fourier transform question This is a question I feel too stupid asking my professor about. I'm having a mental block remembering how this works even though I think I understood it at one point:

I know the following properties:

$$x(t) {\longrightarrow}\boxed{\textrm{LTI System}}{\longrightarrow} y(t) = x(t) \star h(t) \longleftrightarrow X(j\omega)H(j\omega)$$

$$x(t) = e^{(j\omega_0 t)} \overset{\mathcal F}{\longleftrightarrow}X(j\omega) = 2\pi\delta(\omega-\omega_0)$$

So why is this true:

$$x(t) = e^{(j\omega_0t)}{\longrightarrow}\boxed{\textrm{LTI}}{\longrightarrow} y(t) = e^{(j\omega_0t)}H(j\omega)$$

instead of this:

$$e^{(j\omega_0t)}{\longrightarrow}\boxed{\textrm{LTI}}{\longrightarrow} y(t) = 2\pi\delta(\omega-\omega_0)H(j\omega)$$

I know I'm missing something here or have some fundamental misunderstanding, but I can't seem to catch what it is. I'm taking a DSP course, but it's been quite a while since I took basic signals and systems. If anyone could help me out I'd really appreciate it.

You're confusing time and frequency domain. Why should the (time domain) output signal $y(t)$ be a function of $\omega$? If $x(t)=e^{j\omega_0t}$ is the input signal to an LTI system with frequency response $H(j\omega)$, then the output signal is given by

$$y(t)=e^{j\omega_0t}H(j\omega_0)\tag{1}$$

Note that in $(1)$, $H(j\omega_0)$ is a constant, not a function of $\omega$.

Equation $(1)$ can be easily derived from the convolution integral

$$y(t)=\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\tag{2}$$

with $x(t)=e^{j\omega_0t}$ and the Fourier transform relation

$$H(j\omega)=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}dt\tag{3}$$

I'm sure you can fill in the details yourself.

Your last equation would be true for the Fourier transform of the output signal:

$$Y(j\omega)=X(j\omega)H(j\omega)\tag{4}$$

And with $X(j\omega)=2\pi\delta(\omega-\omega_0)$ you obtain

$$Y(j\omega)=2\pi\delta(\omega-\omega_0)H(j\omega)=2\pi\delta(\omega-\omega_0)H(j\omega_0)\tag{5}$$

which is of course just the Fourier transform of Eq. $(1)$.


Determine the system function H(s) of a system and find out the differential equation I have created the following system for practice purposes. From this system I want to determine the system function H(s).

In the picture I have worked with auxiliary (dummy) variables, which should simplify the whole calculation (see blue marking).

How did I approach the matter now? Well, first I have drawn in my auxiliary variables, then I have looked at the summing elements.

I know the relation $H(s)=\frac{Y(s)}{X(s)}$, so I have to transfer my equations into this form.

So these are my equations:

$$\text{Eq 1:}\quad X(s) + G(s) = P(s)$$
$$\text{Eq 2:}\quad 3P(s)=Y(s),$$
$$\text{Eq 3:}\quad G(s) = s^{-1}Y(s)+2Y(s)$$
Eq 3 in 1 results in:
$$X(s) + s^{-1}Y(s)+2Y(s) = P(s)$$
Using Eq. 2 results in:
$$Y(s) = 3X(s) + 3s^{-1}Y(s) + 6Y(s)$$
Therefor:
$$Y(s)(-5-3s^{-1})= 3X(s)$$
$$H(s)=\frac{Y(s)}{X(s)}=\frac{3}{-5-3s^{-1}}$$

My first question is, if you see the system function $H(s)$, is it correct so far?

Then I would be interested in the conversion into a differential equation, here I have already prepared something:

$$Y(s)(-5-3s^{-1})= 3X(s)$$

This should now be handled with the inverse Laplace transformation:

$$-5y(t)+\mathcal{L}^{-1}(-3s^{-1})= 3x(t) \\-5y(t)-3\int_{0}^{\tau} x(\tau) d\tau = 3x(t)$$

My problem at this point is that when I look into my Laplace transform table, I find that $s^{-1}Y(s)$ transformed back results in an integral. Now it is said that these systems can also be given in a differential equation (in terms ot the z-trafo this would be comparable to the difference equation there). How would one come to this differential equation? I have a hunch but I would be interested in your answers :)

I hope my question is so far understandable. Thank you very much.

Your transfer function looks correct. Note that you can rewrite $H(s)$ as

$$H(s)=-\frac{3s}{5s+3}\tag{1}$$

or, equivalently,

$$5sY(s)+3Y(s)=-3sX(s)\tag{2}$$

Since multiplication by $s$ corresponds to differentiation in the time domain you obtain from $(2)$

$$5y'(t)+3y(t)=-3x'(t)\tag{3}$$


Is spectral centroid enough to measure brightness? i'm new in this group and i'm studying a sound signal.
I have used FFT for extrapolate frequencies and see the harmonics of this sound.
I have to compare this sound with another similar sound and find the diffrences.

From various research I found that an important element that characterizes a sound, and therefore it's timbre is the brightness.
My question is:
is it sufficient to calculate the brightness using only the specrtal centroid, using the frequencies and the 'magnitude' of my fft?

'Brightness' is a subjective descriptor of an audio signal. It can't be calculated directly from physical properties of the signal. A panel of subjects can give you their ratings for brightness of one signal relative to another.

That said, the spectral centroid has long been known to be correlated with subjective measures of brightness. So while you are not calculating brightness itself you should be calculating a fairly good indicator.


Why is Fourier space not adequate for (theoretical or digital) filters? As far as I have seen, almost all theoretical filter design occurs in Laplace or Z-space. Also, there is a pervasive connection to real life analog filters in the design. If one is just thinking in a mathematical theoretical thing (or something that could be implemented digitally), why wouldn't one filter signals in Fourier Space?
Why is, say, multiplying the Fourier Transform of a certain function by a unit step up - step down function, and then making the Inverse Transform of the resulted signal a "band pass filter"? Why should one use Butterworth filter or similar things to make a digital filter?

We actually do implement long FIR filters doing that.  If it's running real-time or on an input that is indefinitely long, then the method for doing that is called "fast convolution". Remember that the "Fourier Transform" that you are multiplying is the Discrete Fourier Transform and that has issues regarding circularity in the resulting convolution.
We deal with that problem using "Overlap-Add" or "Overlap-Save".  The frequency response that you are multiplying by must be one of an adequately short FIR to prevent "time aliasing".


Filter, window function and audio signal: how to link them to each other? I have a question about the implementation of digital filters in CoreAudio. I'm in big trouble because it is a few weeks I'm trying to understand how to implement them. The basic idea is this: while I talk to the iPhone's microphone, my voice is filtered by a low-pass filter or high-pass or bandpass.

Reading the book "Digital Signal Processing" and following the algorithms mentioned in the book "Digital Sound & Music" I came to write the following code (high pass filter):

float cutoff = 5000;
float rate = 44100;
int N = 1024;
float f_c = cutoff/rate;
float ω_c = 2*M_PI*f_c;
int middle = N/2;
int i = -N/2;
float fltr[N];

do {

    if (i == 0)  {
        fltr[middle] = 1 - 2*f_c;
    }
    else {
        fltr[i+middle] = -sin(ω_c*i)/(M_PI*i);
    }
    i++;

} while(i != N/2);

After writing the filter, in the book, you are instructed to run a windowing function. Here I have the problem: I do not understand what to do. Reading on Apple's official documentation, I find the function of Blackman.

int N = 1024;

float dataOut[1024];

for (int n=0; n < 1024; ++n)
{
    dataOut[n] = 0.42 - (0.5 * cos(  2 * M_PI * n / N ) ) + (0.08 * cos( 4 * M_PI * n / N) );
}

But once inserted, what should I do? How can I connect the filter carrier, the Blackman function and the audio signal? I just can not understand it.

Thank you for your attention.

UPDATE 1

I answer to users bristow robert-johnson and Matt L. (I write here because the comment was too long to be invited)

Hi, thank you for answering. There is no filtering at this time because I do not understand how to connect the filter (?) With Blackman window.

What I understand from the chapters devoted to the book “Digital Signal Processing” filters: Create a filter, a sort of kernel filter that is a kind of "mask" to be applied to the signal itself. In this state the filter is said "Impulse Response" or also is said IIR. Now I apply the mask to signal by a function convolution becoming a FIR filter. What I have written is this correct?

As for calculating coefficients, as I understand it, they should be the coefficients Nyquist angular frequency. Is right? If yes, should be variables of "f_c" and "ω_c" (this is written in the commentary of the algorithms in the book).

Your sample code looks like it is trying to create a windowed Sinc FIR filter.  A Sinc function (the transform of a "brick-wall" rectangle filter) has infinite support, thus it must be windowed to create a finite length FIR filter kernel.  There is pseudo-code to generate such a filter, using a Hamming window, on my DSP web page.

You can apply such a FIR filter to audio by convolution with your audio samples (either directly convolution, or using FFT/IFFT overlap-add/save fast convolution on blocks of data).

If you are using iOS frameworks such as Core Audio, note that you can also use the vDSP_biquadm API in the iOS Accelerate framework for IIR filtering.


Real-time FIR filtering on time domain My previous question wasn't clear. So I thought it'll be good idea to repost with more details.

I've designed a 100th order windowed highpass FIR filter with 3Hz cutoff frequency for my ADC datas with 100Hz sampling rate. I am currently stacking up my datas in a buffer. Then convolve that buffer with my FIR coefficients array. Lastly, take FFT of the output of convolution. This way works pretty well. I can see clearly the respond that I wanted from FIR filter. However since output is an array, I can't display that to represent my time domain graph.

Since I'm working on a real-time data plotting, I need to send filtered time domain values one by one to show them in my plotting application. To get a single output I've tried to reduce my adc datas' array size to filter's array size and implement $y(n)=h_0x(n)+h_1x(n-1)+\ldots +h_{N-1}x(n-N+1)$ this instead of convolution.

I gave 3Hz 0.2V amplitude sine wave signal with 1V DC offset. I've recorded the result in a video. In the video, you can see it removes the DC offset and amplitude is right but the wave is no longer a sine wave. Later in the video, I've changed the frequency to 1Hz and it cutted off the signal which was expected. In the frequency domain everything is wrong because it shows same frequency values no matter what real frequency is.

The green output in the video looks like it's just the blocks of the filtered data, with transients at the beginning of each block. That's probably why the spectrum always looks the same, you're just seeing the blocking artifact / frequency. You are probably starting the filtering from scratch every time, whereas you need to keep the data (state) of the filter from one block to the next.


How sine wave is generated using lookup tables From this document I am quoting P2,L7

  So what is the incremental distance if we are using a 64 element look-up table instead of a sine function based on 2π? incremental distance = 2k*64/48k = 8/3. In this case increment represents a virtual phase increment or more directly a table index increment.

The radian increment statement I am clear with. Now I am having a basic doubt of how can it be incremented to table index of 8/3 on each sample to be calculated.

Following, some new doubts arise.

Is lookup table method applicable for limited number of Fs, which can be indexed from the lookup table ?
Consider a lookup table of size 256, generated with resolution of π/24. Is it popssible to generate sine value for, say sin(π/3) from lookup table ?

Thanks in advance

Lookup tables are incredibly useful.

In audio, when you're doing a chorus or flange effect, or a synthesizer with a modulating envelope, you want to be able to cycle through a sine wave in a lookup table and not just keep track of some infinitely growing index n.

You'll want to use some mathematical program like MATLAB, Octave, or Mathematica to generate as many points as is practical (I wouldn't go with less than 512 or 1024 points) and store that in an array. This is your lookup table.

For a fractional index, what you'll end up doing is linearly interpolating around the nearest two indices at that fractional value. Let's say your index is 8/3. You've called up your index in a 1024 point table four times. What do you do? Now you're at index 32/3.

You'll need to find the lower index floor(32/3) and the higher one ceil(32/3). You'll need some conditionals in case you're at the end. I'll leave that exercise to you.

So now you've got your neighbor indices, you'll evaluate the interpolated value at that fractional index using the linear interpolation formula.

y = x*y2 + (1-x)*y1;

where x is the difference between your fractional index and the nearest index rounded down, y1 is the nearest value of the lookup table rounded down (you found the index earlier, now just use it to index into the lookup array) and y2 is the same thing for the nearest value rounded up.

You would track the index in a variable and increment it each loop of your program by 8/3 and do a wraparound when necessary (when you go beyond the length of the lookup table).

The trickiest thing is keeping track of the indices. The nearest value rounded up may be the 0th index.


Signal Resampling by a non integer factor with convolution to start I am a newbie in signal processing, I just started a month ago so please be as detailed as possible.

From what I understood to resample a signal by a non integer factor you can do an upsample followed by a downsample. So, imagining a sample set of size 10 and if you want to upsample it to 15, you can upsample it by 3 and downsample it by 2, correct?

If that is true, does one can create a kernel by merging the two kernels the upsampling and downsampling ones?

To upsample by 3 we can use the kernel [1/3, 2/3, 1, 2/3, 1/3] and to downsample by 2 we can use the kernel [1, 0]. Can I create just a kernel to do both operations?

If instead I use interpolation, can I do it using a kernel as the above?

You might be describing something similar to a canonical polyphase resampling filter.  If you can compute your filter+interpolation kernel on the fly (certain windowed Sinc’s and others) you can even do this for irrational sample rate ratios or factors.


When does a Window function's sidelobe level become out of band attenuation? Another question I'm struggling to answer by myself.

For a Window function, say Blackman (response pictured), what is the out of band attenuation? If I'm designing a filter with out of band attenuation of 40dB, is it the level of the first sidelobe that I'm interested in? Or is the out of band attenuation dictated by something else (like 3rd sidelobe level for example?)

Thanks!

Windows with low sidelobe levels have large transition
bandwidth's. Transition bandwidth is inversely proportional to N for a
given window.  Indeed, the ratio of transition width over sidelobe level
in dB for fixed length N is virtually a constant.

We can use the following empirical formula to estimate the required filter length

where,

N is filter length,
A is the out of band attenuation in dB
fs is the sampling frequency


Accuracy of Goertzel algorithm for measure harmonic If I use Goertzel algorithm for calculate the magnitude of current or voltage, use for measure current harmonic or voltage harmonic or fundamental component. I not use for detect harmonic. Do you know the accuracy for measure harmonics is good or not?,How much about relative error?

Another one question is The Accuracy relative with number of data sample?
Ex. If I have measure voltage harmonic,
 I measure waveform of current by use sample rate = 4000 Hz, target frequency for measure fundamental harmonic is 50Hz, so if I use number of sample (bin width) for Goertzel algorithm is 80 sample this is minimum, If I use more then 80 samples like 4000 sample the accuracy will more accurate right? how much better?

The most common Goertzel algorithms are identical to calculating the magnitude of 1 bin of an FFT/DFT using a rectangular window.  As such, there will be rectangular window artifacts (so called "leakage") for any spectral content, such as harmonics, that are not exactly periodic in the DFT length, or Goertzel length.

Thus the accuracy of using just 1 bin may not be more accurate if you increase the length from one where all the spectral content is periodic to one where it isn't.  With a complete FFT, one can somewhat mitigate this effect by using high quality interpolation in the complex domain before taking the magnitude.


Should we use a digital low pass or high pass filter to remove sensor bias? Consider a continuous signal oversampled at, say $2 \;kHz$, and then system digital low pass filtered to a $100\;Hz$ frequency which is the control loop frequency. It is known that there is some bias in the signal since it is an accelerometer output.

If we use a discrete derivative with a discrete integral in series, will this remove (atleast some of) the DC bias in the signal?
Correspondingly it seems a better idea to use a high pass filter rather than a low pass filter although this will invite noise. Is a notch filter how people deal with accelerometer signals and should one go with a high pass approach in case of the unavailability of a notch filter (since the signal necessarily contains bias)?
Does the same logic hold if one designs the filter in the continuous domain?

If we use a discrete derivative with a discrete integral in series, will this remove (atleast some of) the DC bias in the signal?

Yes (if your integrator starts at state 0).
Note that this is a CIC filter.
Also notice an offset should be present in your accelerator data, or you just got lost in space.

Correspondingly it seems a better idea to use a high pass filter rather than a low pass filter although this will invite noise.

Well, seeing that by definition of "low-pass filter", the DC component remains untouched, this doesn't only seem to be a good idea, but mandatory.

Is a notch filter how people deal with accelerometer signals and should one go with a high pass approach in case of the unavailability of a notch filter (since the signal necessarily contains bias)?

It depends on what people are doing for what purpose. My guess is you want to estimate something from your accelerometer observation. Maybe Kalman filters are what you're after?

Does the same logic hold if one designs the filter in the continuous domain?

This has nothing to do with continuous or discrete time. So, yes.


Problem designing a specific filter I have the next problem.

$H_{c1}(j\omega )$ is the ideal antialising filter and $H_{c2}(j\omega )$ is a real one. I'm asked to design $H(e^{j\Omega })$ so that $y[n]$ in the second diagram (the one on the right) is exactly the same as $x[n]$ in the first diagram (the one on the left).

I did something but I'm pretty sure that it isn't right. I found the equation for $H_{c2}(j\omega )$ as a function of $\omega$. Then I normalized the $\omega$-axis due to the A/D converter and divided the magnitude by $T$. If I didn't make any mistakes while doing the algebra, we get that

$$H_{c2}(e^{j\Omega })=\frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \Omega + \frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \pi + \frac{1}{10}\hspace{0.5cm} for\hspace{0.5cm} -\pi < \Omega < -\omega _{p}T$$

After that, I just thought of finding the inverse function of $H_{c2}(e^{j\Omega })$ for $-\pi < \Omega < -\omega _{p}T$... And that was it. Let $B(\Omega)=\frac{1}{\frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \Omega + \frac{9}{10}\cdot \frac{1}{\pi - \omega _{p}T}\cdot \pi + \frac{1}{10}}$. Then,

$$H(e^{j\Omega }) =
\left\{
 \begin{array}{l1}
  B(\Omega )  & \mbox{if } -\pi <\Omega <-\omega _{p}T \\
  1 & \mbox{if } -\omega _{p}T < \Omega<\omega _{p}T \\
B(-\Omega) & \mbox{if } \omega _{p}T < \Omega<\pi
 \end{array}
\right.$$

The problem with this is that I'm pretty sure that the filter I "designed" has no antitransform. I think that there is a smarter way of approaching this exercise but I just can't see it.

so that $y[n]$ in the second diagram (the one on the right) is exactly the same as $x[n]$ in the first diagram (the one on the left).

Considering that $H_{c1}$ is per your very own definition $\ne H_{c2}$ means that without restricting your system's input $r_c$, there's nothing you can do. You will have to implement a filter with transfer function $H_{c1}$ if you need a system with transfer function  $H_{c1}$; there's no way around that.

Now,  $H_{c1}$ being impossible with finitely long filters (as it is infinitely sharp in frequency domain), this is impossible.