Signals Cross Correlation from Their Power Spectral Densities? For two real valued signals $f$ and $g$ with their corresponding power spectral densities $|F(k)|^2$ and $|G(k)|^2$ Is there a way to get an intuition of how correlated $f$ and $g$ are just by looking at the power spectral densities (PSDs)?

For example does the fact that the PSDs of these two signals are so similar imply that they are cross correlated?

So far I have considered that $$|F(k)|^2 = F(k) \overline{F(k)}$$  $$|G(k)|^2 = G(k) \overline{G(k)}$$

and the cross correlation we want to find is

$$ F[f \star g] = F(k)\overline{G(k)} = G(-k)\overline{F(-k)}$$

and so it seems that there could be a way to get an exact result for the cross correlation from the PSDs, however I would also be happy just to get a hand wavy  idea of if these two signals are cross correlated

No,

For instance 2 instances of AWGN with the same STD have the same PSD while being totally uncorrelated (Moreover, Statistically Independent).

Also, pay attention that Correlation is not commutative:

$$ f \ast g \neq g \ast f \Leftrightarrow F \left( k \right) \overline{G \left( k \right)} \neq G \left( k \right) \overline{F \left( k \right)} $$


Convolution from bottom right I want to do a convolution from the bottom right and not as usual from the top left. I think conv2 of Matlab only does from the top left.

How can I do a convolution in Matlab from the bottom right?

Thank you very much for the answers.

I hope I understand your question correctly, in that you're trying to produce a mirror image of the convolution kernel (filter) and then convolve. Flip you convolution kernel. In MATLAB, you can use the flip command. If you flip it left-to-right, this should do it.

However, if you're simply saying that you want to use the same kernel but start the convolution process from a different point, this would not give you a different result back at all. This is one of the reasons you can do FFT based convolution.


Wavelet thresholding What is the difference between soft thresholding and hard thresholding. Where we use soft and hard thresholding in image for denoising. I understand that in hard thresholding, the coefficients below threshold value are set to zero and the value above the threshold is set to one. Please explain me about soft threshold.
Please explain whether the threshold value is the intensity value of the image. For example if the intensity value ranges between 0 to 255. In case of hard thresholding if the threshold value is considered as 100 then the values below 100 is set to 0.The value above 100 are retained. Is this correct? Please correct me.

For a given threshold $\lambda$ (that can be dependent on resolution level), and value of wavelet coefficient $d$, hard thresholding is defined as:

$  D^H(d|\lambda)=\begin{cases}
               0,& \text{for } |d| \leq \lambda\\
               d,& \text{for } |d| > \lambda
            \end{cases}$

whereas soft thresholding is governed by following equation:

$  D^S(d|\lambda)=\begin{cases}
               0,& \text{for } |d| \leq \lambda\\
               d-\lambda,& \text{for } d > \lambda \\
               d+\lambda,& \text{for } d < -\lambda \\
            \end{cases}$

Figure below depicts both cases:

The soft thresholding is also called wavelet shrinkage, as values for both positive and negative coefficients are being "shrinked" towards zero, in contrary to hard thresholding which either keeps or removes values of coefficients.

In case of image de-noising, you are not working strictly on "intensity values", but wavelet coefficients. You probably remember that you can decompose your image into wavelet levels, like in case of lovely Lena.

 Assuming that wavelet transform gives sparse coefficients, mostly close to zero, and noise level is lower than wavelet coefficients, you can simply threshold these. Although if you wish, you can perform hard/soft thresholding on each decomposition level with a different value of $\lambda$. When it is done, then you just have to reconstruct your image from all decomposition levels and voila, noise should be removed!

Below you have two examples of de-noised image via hard and soft thresholding respectively (same $\lambda$). Obviously soft thresholding gives more smooth image- if you can notice that with such a poor resolution ;) Courtesy of MATLAB.


Linearity with Difference Period in Fourier Transform I know that a system is linear if it satisfies

$$\mathscr{F}\{ a\,x(t)+b\,y(t) \} = a\,X(\omega)+b\,Y(\omega)$$

for Fourier transform, $X(\omega)\triangleq\mathscr{F}\{x(t)\}$

But what if $x(t)$ and $y(t)$ don't have the same period?

Can we use same linearity condition in Fourier transform?

Continuous time Fourier transform does not require the signals to be periodic, and therefore the linearity assumption holds. CTFT is a linear operator.

Furthermore, if you mean CT Fourier series, which require that both signals to be periodic by (a common period) $T_0$, that's also linear.

But if input signals have no common periods, then you cannot apply the linearity test for CTFS.


Detecting phase difference between a clean and a noisy signal I have a clean square-wave signal and a noisy signal from an accelerometer,
perfectly synchronized to each-other but out of phase. I want to find the phase difference between them in real time, in software (I am building a dynamic balancer).
I plan to convert the clean signal to a (complex) phasor in real time and multiply that by the noisy signal, then integrate over a moving window. I should be able to measure the phase angle of the integral to get the answer.
Does this sound like the right approach? It would seem to cancel out the other frequencies automatically.
I have considered measuring the frequency of the clean signal and filtering the noisy signal using a band-pass filter around that frequency, but it seems unnecessary.

Here is how I would do it:

1) Find the frequency of the clean signal, this should be trivial in the time domain for a clean square wave.

2) Heavily smooth the noisy signal with a symmetric filter.  (IIR or FIR)

3) Calculate a single DFT bin for each signal using the length of a whole number of cycles for your frame.  This will make the harmonics leakless, and make your fundamental measurement more accurate.  The bin number you calculate should be the number of periods in your frame.

4) Calculate the phase for each signal relative to the frame at the fundamental frequency.  (Straight from the DFT bin value)

5) Take the difference in the two phases and translate if from an angle to a time using the frequency.

This is going to be way more efficient than some cross-correlation type calculation.


Using an error prediction filter for filtering a semi-known signal I'm trying to wrap my head around the proper use of a Wiener or error-prediction filter for filtering data. It seems to me that it is only a whitening filter, so how is it used when the data you want to recover isn't an AWGN signal?

For instance, I have a signal that has several distint interfering signals - I can see them on a PSD, but I don't know that they are a) stationary and b) what properties they have. I can use a method like Yule-Walker equations to recover the AR model for the whole signal, but in this case I only want to recover the model of the interfering signals, not the portion I want to recover.

I tried implementing an adaptive LMS notch filter, with the reference signal being a single sinewave, but this turned out to me much too narrow and didn't track frequency changes in the signal very well.

I guess basically my question is this, if I'm using an error prediction filter to filter real data, then how do I separate the data portion from the noise portion? In other words, I don't want to whiten the whole signal, only the noise portion. What am I missing?

i am not sure i understand correctly the question (feel free to update me if not so).

There is the MUSIC algorithm, which extracts signals embedded in a background noise, as a sum of sinusoidals signals

There is also the option to use SVD (or Karhunen-Loeve transform) and reduce the dimensionality of the input data while retaning maximum information (this will discard most of the background noise components).

If this is online or real-time, this could be done adaptively.

Hope this helps


Performing classification based on FFT results It is perhaps important to start with the fact that I am a complete beginner in DSP.
I have got a number of audio recordings (of sampling rate 22 kHz) - of bird songs - which I have been trying to analyse using FFT (in Matlab/Octave). Particularly, I am trying to show using Machine Learning classification algorithms that different classes of those recordings have prominence of different frequencies (or frequency ranges).
The recordings are of variable length and, due to computation limitations, the largest size FFT that I can do is 2^19 (which I understand is the number of points it takes from each audio file).
So, my first question is: if I break my recordings in parts, each corresponding to the size of FFT that I've chosen - would it be still reasonable to treat those parts as separate data examples (i.e. separate recordings), and what kind of information do I lose when splitting the larger recordings in such way?

The second question is: is there a better way for a beginner to perform this analysis, in a not so computationally expensive way, since I think working with vectors of size 2^18+1 is not really the best thing to do in the current case.

I think the FFT is a bad choice of representation for your problem - it captures many properties of the signal irrelevant to your application, and as you are suspecting, it generates a huge amount of data to process if you extract the FFT of the whole signal.

It seems to me that the most important quantity to consider when studying birdsongs is pitch (fundamental frequency) - all other dimensions of sound (loudness, timbre) are actually variability factors you want to get rid off. For example, two recordings of the same birdsong made in a different environment and with different equipment will exhibit a different frequency response due to the variation in conditions; but fortunately the pitch profile would be exactly the same!

So I suggest you to use a pitch transcription utility (canned solutions : aubio, praat, sonic visualizer...) to extract a pitch contour - a function giving the predominant "note" as a function of time. From that, you could define a feature vector containing pitch statistics (mean, standard deviation, maybe higher order moments); or maybe just build a histogram of pitch values; and this would result in a very compact feature vector suitable for automatic classification. To improve your results, you might then add features capturing the dynamics of pitch over time - dominant modulation rates extracted form the pitch contour, variation of pitch statistics over chunks of a few seconds of audio, etc..


How to choose the size of a Laplacian of Gaussian kernel for filtering images fast? I am filtering an image using a Laplacian of Gaussian (LoG) kernel. The kernel dimensions are the same as the image. The filter operation takes very long to complete. How can I speed this up while maintaining fidelity?

The LoG kernel is not separable like the Gaussian kernel. So I cannot speed it up by filtering the image with 2 vectors successively. The only other option is to resize the kernel to get rid of insignificant multipliers on the outside. Is there a rule of thumb when choosing the size of an LoG kernel (considering image size and standard deviation of LoG)?

I am using this picture (633x900):

Filtering this in matlab using Gaussian vs. LoG gives vastly different times with the same kernel size:

>> I = imread('PICTURE.jpg');

>> tic;imfilter(double(I), fspecial('gaussian', [633,900], 50)); toc
Elapsed time is 1.129423 seconds.

>> tic;imfilter(double(I), fspecial('log', [633,900], 50)); toc
Elapsed time is 253.447840 seconds.

I just started with image processing so please excuse me if I missed out on something glaringly obvious. Thanks!

you can approximate LoG with Difference of Gaussians. Apply 2 Gaussians (separable) with different sigmas. (for what I know with sigma ratio 1.1) and take the difference.


Estimate phase using STFT in Matlab? My signal oscillates back and forth on the interval [-1, 1] at varying frequencies, with some added higher frequency noise components (example data shown below). I want to to estimate the phase of the slow oscillation at any given time. This means I would like to filter out some of the high frequency noise, and then calculate the phase at each time for the smoothed signal. However, my method isn't working so I'd really appreciate any ideas about how to do this properly!

I've attempted to do this with the Short Time Fourier Transform (FTST). The method I've tried is to calculate the Fourier transform for a sliding window across my signal, choose the frequency component with the highest magnitude, and then use the phase from that component.

However, I'm having trouble getting the result to match up with the slowly oscillating component that I can see in the signal. I've tried two methods and neither looks good:

using just the phase shift from the complex value of the fourier transform computed using Matlab's spectrogram() function
using the total phase based on the time and frequency return variables from the STFT plus the phase shift

Am I missing something? I'm new to signal processing, so it might be a conceptual misunderstanding or just an error in my implementation. I would really appreciate any help or suggestions!

Here is the Matlab code I'm using to generate the figure:

% define STFT options
window = 50;
num_overlap = window - 1; % because we need phase at every position
signal_length = length(signal);

% pad signal with zeros on both sides
pad_left = floor(num_overlap/2);
pad_right = num_overlap - pad_left;
padded_signal = [zeros(1,pad_left) signal zeros(1,pad_right)];

% calculate STFT
[S,F,T] = spectrogram(padded_signal,window,num_overlap);
component_magnitude = abs(S);
[m max_index] = max(component_magnitude); % get the frequency with maximal amplitude

% calculate the phase when the frequency components are maximized
for ti = 1:size(S,2)
   phase_shift(ti) = phase(S(max_index(ti), ti));
   signal_phase_total(ti) =  T(max_index(ti))*F(max_index(ti))+phase_shift(ti);
end

num_plots = 5;
% plot the initial signal
subplot(num_plots,1,1)
plot(signal, 'LineWidth', 2);
title('signal', 'FontSize', 18)
axis tight

% plot the STFT
subplot(num_plots,1,2)
imagesc(component_magnitude);
hold on
plot(1:length(max_index), max_index, 'y', 'LineWidth', 2)
hold off
colormap(bone);
title('STFT magnitudes with max component in yellow', 'FontSize', 18)

% plot the phase info
subplot(num_plots,1,3)
plot(phase_shift, 'k');
title('phase shift of maximal component', 'FontSize', 18)
axis tight

% try to superimpose the phase shift onto original signal
subplot(num_plots,1,4)
hs = plot(signal, 'LineWidth', 2);
hold on
% plot the phase of the max component
hp = plot(real(cos(phase_shift)), 'r');
lh = legend([hs hp], 'Signal', 'Phase');
set(lh, 'FontSize', 12);
axis tight
hold off
title('phase shift of max component superimposed onto signal', 'FontSize', 18)

% try to superimpose the total phase onto original signal
subplot(num_plots,1,5)
hs = plot(signal, 'LineWidth', 2);
hold on
% plot the total phase
hp = plot(real(cos(signal_phase_total)), 'r');
lh = legend([hs hp], 'Signal', 'Phase');
set(lh, 'FontSize', 12);
axis tight
hold off
title('total phase superimposed onto signal', 'FontSize', 18)

Jumping between bins of STFTs will lead to the discontinuities seen in your phase estimates.  Interpolating both frequency and phase from a sequence of STFTs is possible, but requires non-standard techniques, such as fftshifts and 2D interpolators.

If you know the frequency range of your desired slow oscillation, then a bandpass filter in conjunction with a PLL follower might provide a better looking phase result using more standard DSP techniques.


How to find the instantanous SNR over Fading channel I would like to prove that the instantnous SNR of flat fading channel
$$y=hx+n$$
is
$$\gamma=\frac{|h|^2}{N_0}P$$

where $E\{z\}$ is the expectation of $z$ and

$$
E\{x\}=P
$$
$$
E\{n\}=0
$$
$$
E\{n^2\}=N_0
$$

We know that $h$ and $n$ are random varibles so

The noise power is
\begin{align}
E\{n^2\}&=N_0
\end{align}

The problem for me is in the Received signal power

\begin{align}
E\{|hx|^2\}&=E\{hxh^*x^*\}\\
&=E\{|h|^2\}P
\end{align}

Why they take $E\{|h|^2\}=|h|^2$ and we know that $|h|^2$ is random varible and it has it own expectation.

Also what if a noise $n_1$ is given by

$$ n_1=hn+n$$

what is the noise power of $n_1$ where h is the fading and $n$ is nose with zero mean and power $N_0$.

Thanks.

You are trying to find the average SNR. In instantaneous SNR, the power of the channel is $|h|^2$.

Your SNR equation isn't accurate. $N_0$ is usually denote the PSD of the noise. The noise power is $\sigma_n^2 = N_0\,W$, where $W$ is the signal bandwidth. Or you define SNR as $E_b/N_0$ where $E_b$ is the signal energy.


DTFT Pairs confusion When I am in the DT Fourier Domain, and I want to come back to the time domain, which pair do I use? Asking because both pairs have the exact same "form" in the Fourier domain, and that is what is tripping me up. Do I just use pair 2 when I want to go into the Fourier domain and then I have to use pair 1 when I want to come back?? Or is there more to it than that?

I would really recommend that you compute the DTFT of the sequences $x_1[n]=a^nu[n]$ and $x_2[n]=-a^nu[-n-1]$. This is very straightforward if you use the formula for the geometric series:
$$\sum_{n=0}^{\infty}q^n=\frac{1}{1-q},\qquad |q|<1\tag{1}$$
It will be important to consider the magnitude of $a$. Note that it's not for nothing that the table in your question specifies a condition on $|a|$.


`SSIM(YUV)` vs. `SSIM(Y)`: Values Are Very Different I'm seeing some weird results when trying to compare image qualities in a few different color spaces. Have a look at my graph:

I want to compare with another publication, which reported PSNR(YUV) and SSIM(Y). However, I don't understand how I can have SSIM(YUV)=0.9 for the highest quality point on the right, but only have SSIM(Y)=0.7 for the exact same image. How can there be such a big difference? Is this normal? Or am I doing something wrong?

For completeness, here are reference image and the image under test:

And the MATLAB code I use for converting RGB to YUV:

function YUV=rgb2yuv(RGB)
R = RGB(:,:,1);
G = RGB(:,:,2);
B = RGB(:,:,3);

Y = 0.299   * R + 0.587   * G + 0.114 * B;
U = 0.5 - 0.168736 * R - 0.331264 * G + 0.5 * B;
V = 0.5 + 0.5 * R - 0.418688 * G - 0.081312 * B;
%   ^^^ I tried to use 0.0 as well, but SSIM values are identical

YUV=cat(3,Y,U,V);
end

I figured it out: the ssim() function in MATLAB, does not handle RGB images, but instead threats it as a 3D-image (with the third dimension time for example, to measure the quality of a video). So it interpreted my 3-color-component image as a video with 3 frames.


Masking the mains hum to improve anonymity Today, a news blurb appeared that claimed that attackers could identify the location in which an audio record was taken by matching the mains hum to known data. Sadly, I fear that this might be feasible. The mains hum has quite a tight tolerance, but any aberrations (that cannot be done away with) could be spotted easily. In fact, power suppliers have to observe the hum and adjust their stations accordingly.

How could you manipulate a record such that this information is destroyed completely, so that it is not reconstructible by any means? Filtering alone might leave enough remains to still perform that identification. After all, all filters have a finite number of poles.

I believe it's possible to get location information from 50/60 Hz hum, particularly if the recording is long-duration.  See this paper by Catalin Grigoras from 2009.

One way to defeat this would be to add enough sampling rate jitter to your recording so as to swamp the small variations in mains hum.  That paper showed information from variations of around 1 part in 1000, with sampling every 10 sec or so.  If you let the sampling rate of your signal drift around randomly on that kind of scale (which would have no perceptible effect), it would obscure the variations in the hum.

So, for instance, if you let your sampling rate drift by +/- 0.1% (e.g., 44,050 to 44,150), with a maximum rate of change of, say, 1 Hz in 10 sec, that would impose a random component on top of the underlying power-grid-derived variations, making them hard to separate.  This modification would be completely inaudible.


How does causality (i.e. unit step) affect the DTFT of a sine or cosine wave? Tables of common Discrete-Time Fourier Transform pairs list the transform of a sine wave:

$ \sin(\omega_0\ n) $ and its transform:
$ -j\pi\ [d( \omega\ - \omega_0\ ) - d( \omega\ + \omega_0\ )] $

And the cosine:

$ \cos(\omega_0\ n) $ and its transform: $ \pi\ [d( \omega\ - \omega_0\ ) + d( \omega\ + \omega_0\ )] $

How might the results differ if the sin or cosine is causal? That is, I would like to determine the Fourier Transform of the following (sin or cosine multiplied by the unit step):

$ \sin(\omega_0\ n) u[n] $ or  $ \cos(\omega_0\ n) u[n] $

So regarding the DTFT of the signals above, is the only method to realize that multiplication in the time domain is convolution in the frequency domain? Or is there a more simple rule-of-thumb, or property, of the Fourier Transform that I am missing?

(This is not homework, but it is for studying, so it would be helpful to include resources along with answers.)

Multiplication in the time domain is convolution in the frequency domain. The frequency response of a unit step has 1/frequency shape and it's pretty messy.

What you will see is that the overall frequency response is not a delta impulse anymore but that the lines widens and that the get "skirts". These will bigger and more "wiggly" for the cosine than for the sine as the step function slices right through the maximum of the cosine.


What are True Positive and False Positive with reference to image processing Say for instance, face detection algorithms - I see researchers testing their algorithms and table their results under TP and FP. What are these?

A ,,true'' positive would be an image that has a property (in the context of face detection maybe: ,,contains a face'') and that is recognized by a program as such. A ,,false'' positive then is one that does not have the property but is recognized anyway.


Resampling with and without replacement for estimating significance of spectral components In order to test the significance of spectral components, it seems reasonable to randomly sort the data in order to destroy all the serial correlations / spectral order e.g. 100,000 times, and then estimate the significance or C.I.s by comparing your sample spectrum against the spectra of all the resorted data. Essentially H0 = there is no serial correlation / spectral order, and therefore the signal is white noise.

Here the resorted data would be resampling without replacement. My question is - is this correct or incorrect for spectral estimation? Would it be better to use resampling WITH replacement? And, if so, then why?

Thanks all for your responses! Indeed I am aware it is called bootstrapping - and just to clarify, my intention is not only to resort the time series 10,000 times, but to also create 10,000 spectra from these resorted versions of the time series in order to establish how many times my spectral component of interest is above the resorted time series spectra.

So H0 is that, over 10,000 comparisons, my spectral component is not higher in amplitude than the corresponding component in the resorted spectra more than 5% of the time. My H1 being that the spectral component is greater in amplitude (and therefore not due to chance) in at least 5% or more of the comparisons.

Indeed, I think you can go further than this and also calculate how many times it is BELOW the resorted time series spectra - this would, presumably, give you an indication of its persistence (i.e. it is more persistent than white noise and is therefore unlikely to change in a given period).

The code below calculates: the spectrum for a sine wave X[k] (length 30 time points), an average spectrum for the sine wave resampled 20000 times without replacement, an average spectrum for the sine wave resampled 20000 times with replacement. The original signal spectrum could as well be an average spectrum of, say, 200 signals of interest.

 Fs=1; %Sampling Frequency
 T=30; %Total signal duration in seconds
 t=0:(1./Fs):(T-(1./Fs)); %Time vector
 f1=1; % A frequency figure (in Hz)
 x=cos(f1*t);
 X=fft(x);
 X = abs(X);

 for i = 1:20000
     ps = round(rand(1,30)*29)+1;
     x_p_all(i,:)=(x(ps))';
     X_hat(i,:) = abs(fft(x_p_all(i,:)));
 end

 X_hat_mean = mean(X_hat);

 for i = 1: 20000
     ps2 = randperm(length(x));
     x_p_all2(i,:)=(x(ps2))';
     X_hat2(i,:) = abs(fft(x_p_all2(i,:)));;
 end

 X_hat_mean2 = mean(X_hat2)

 plot(X_hat_mean)
 hold on
 plot(X)
 plot(X_hat_mean2)

It can be seen that the resampling with replacement averages to be the same as the resampling without replacement average (as expected). Obviously, the energy of an individual resampled signal will be different for with and without replacement, but bootstrapping undertakes the process thousands of times.
Now, if we don't permute and take the average, we have the option to compare our spectrum against each individual permutation, where each comparison would be a comparison like this:

Clearly, in this instance, the spectrum of the resampled signal (in this case, without replacement) looks far from white noise because it contains spurious periodicities as a result of the short time duration and limited number of frequency bins. My method compares these resampled signal spectra against the original signal spectrum e.g. 10,000 times to establish significance - a similar method to bootstrapping (which does not, as far as I am aware, require Bonferroni correction).

...is this correct or incorrect for spectral estimation?

From what is provided in this description, I am more inclined to say that it is incorrect.

Here is my understanding of what you are describing:

You have some discrete time series $x[n]$ and its spectrum $X[k]$ which I will assume here that you obtained via the Discrete Fourier Transform (DFT) as $X = \mathcal{F}(x)$. $n$ denotes discrete time, $k$ denotes discrete frequency. Now, you pick some spectral component $X[k]$ and ask yourself: "Could this component have emerged by chance or is it really part of the signal?".

To answer this question, you shuffle $x[n]$ through a series of permutations that swap around some $n+u$ sample to $n+v$ sample. After shuffling a sufficiently long time, what you end up with is some $\hat{x}(n)$, accompanied by some $\hat{X}(k)$, which is essentially a white noise sequence using the same values of $x(n)$.

And then what? Where are the data for your $H_0$?

I am assuming that you are interested on the significance of some $X[k]$ (?).

If you create some distribution of spectral components from $\hat{X}(k)$, you could answer the question "Does $X[k]$ come from this distribution?". That is, could it be because of chance?

BUT, that would be equivalent to considering the average amplitude of a white noise sequence with equal power as your $x[n]$. Then, if $X[k]$ is above that average you can be sure that it cannot be due to noise.

Ideally, what you would like to have had is a distribution for each X[k]. That is, many realisations of the process that generates your signal, many realisations of their spectra and from all of these, one distribution for each $k$. Then, you could answer the question "Given the value of this $X[k]$, does it really come from my signal or not?"

If you are worried about the variance of the components, you might want to use Welch's Method.

  Would it be better to use resampling WITH replacement?

That would modify the permutation rules to allow or disallow $u=v$. I don't see where is the re-sampling here. You had one spectrum before, you have one spectrum after the re-shuffling...Except if you mean that you would be storing each and every permutation result (?) in which case, yes, that would be wrong because it would represent the statistics of the deterioration process....Or maybe I am getting the whole thing wrong (?).

By the way, there is another thing you need to be careful of. If you did a number of statistical tests to test each and every spectral component then you would be assuming that spectral components are independent of each other and also you would have to correct your scores for repeated tests.

Hope this helps.

EDIT:

A few notes with illustrations to understand better the problem and how to deal with it:

Let's build a signal to work with:

Fs=8000; %Sampling Frequency
T=4; %Total signal duration in seconds
t=0:(1./Fs):(T-(1./Fs)); %Time vector
p=2.*pi.*t; % Phase vector (These last two vectors save us typing the same things again and again).
f1=440; % A frequency figure (in Hz)
%
%
x=cos(f1.*p); %The signal
ps=randperm(length(x)); %A random permutation
x_p=x(ps); %The permuted signal
X = fft(x); % The spectrum of the signal
X_hat = fft(x(ps)); %The spectrum of the permuted signal

The key "objects of interest" here are x (the signal) and ps which is a permutation sequence. The permutation sequence will "loop" through the samples of x and re-map them to different time instances.

Randperm samples WITHOUT replacement. If you sample WITH replacement, then the power of the permuted sequence will be different to the power of the non-permuted sequence (see below). Even if you pass x through ps once, you will get a totally randomised x but you can try doing it a few times anyway.

Here are spectra of x and x[ps[n]] (i.e. x permuted once):

Here is the code that produces that graph:

subplot(221);plot(t,x);xlabel("Discrete time (s)");title("x[n]");grid;
subplot(222);plot(t,x_p);xlabel("Discrete Time (s)");title("Permuted x (x[ps[n]])");grid;
subplot(223);semilogy(Fs.*(0:length(x)-1)./length(x),(X.*conj(X))./(length(x).^2));grid;title("X[k]");xlabel("Discrete Frequency (Hz)");
subplot(224);semilogy(Fs.*(0:length(x)-1)./length(x),(X_hat.*conj(X_hat))./(length(x).^2));xlabel("Discrete Frequency (Hz)");title("Spectrum of permuted x[n]");grid

Notice the normalisation by length(x). Also notice how x is permuted once by x(ps) which effects the "re-mapping".

Is the power of x different than the power of x(ps). No, but only if you sample without replacement, in which case, each instantaneous value of x is re-used exactly once.

If you sample with replacement, some values of x might be re-used (re-map the same time instance twice), in which case the sum will work out to be different. You can try it with some other permutation sequence like ps = round(rand(1,length(x)).*length(x)); which as you can see uses rand which doesn't care about its "history" of random values produced.

Here is the first thing that I am trying to say / clarify from what you are writing:

If you pick one X[k] and compare it against $\hat{X}$ which is the whole spectrum of the permuted sequence (fig above), then:

What you would really be asking would be "Can this component have risen by chance from a random sequence with equal power as my x?". Irrespectively of whether that bottom right spectrum is the result of one permutation or more it will be the spectrum of white noise of a sequence of power equal to that of x. Then, why go to the trouble of permuting x a few times? What is the added value? What is the power of x? The power of x is sum(x.^2). What is the average power? The average power is sum(x.^2)./length(x). How many points are there in the FFT? There are length(x) points. What is the sort of power we expect each bin to take? (sum(x.^2)./length(x))/length(x) = sum(x.^2)./(length(x).^2). This is what the bottom two graphs depict (adjusted for one / two sided spectra, so, basically, each line divided by two).

This is why I am writing above, if you want to proceed this way, then you can simply say, the closer my spectral line gets to that "theoretical" average white noise value, the less I can be sure that it is due to signal.

You mention confidence intervals. If you proceed with the above frame of mind, you would be putting confidence intervals on spectral lines obtained by a distribution across the spectrum. In other words, you would put confidence intervals on $\pm X[k]$ from $\hat{X}$ (?). If you want to put confidence intervals on $X[k]$, you need to look within $X[k]$. In other words, put confidence intervals on the variation of $X[k]$ itself. How can you do this? Instead of deriving your $X$ from $x$, you derive it from successive overlapping frames obtained by $x$. In this way, you form $X_0[k], X_1[k], X_2[k], X_m[k]$ where $m$ is the $m^{th}$ frame and $k$ is the discrete frequency component. So, now, combined with the assumption that $x$ is Gaussian and that the power is squaring a Gaussian therefore the distribution will be Chi squared (the DFT is linear so no further changes), you can put confidence intervals from a Chi-square distribution with the degrees of freedom equal to the number of frames minus one. (Which, by the way, is the way pwelch does it and why I am hinting at Welch's method above).

That is all I am trying to say.

With 30 time instances, you might be able to average "across subjects" (?). So, if your study has 30 participants, you can derive the "normal" variation of each $X[k]$ from those. Then, for an unknown spectrum, you can use those distributions to judge if a specific harmonic is "expectable" there or not.

Hope this helps.

EDIT 2:

The above code can be used to experiment further, especially with re-shuffling. In Octave (and possibly MATLAB too), you can apply the reshuffling more than once in two ways, either via x(ps)(ps)(ps)(ps)(ps) and so on, or x(ps(ps(ps(ps(ps(ps(ps))))))) which permutes the permutation of the permutation of the permutation [...] of x. I have a sneaking suspicion that the first way doesn't work in MATLAB (that is the x(ps)(ps)...(ps). Of course, you can even do it with different permutation sequences every time too.


Why are the observation features of an HMM-based recognition/synthesis system modeled by a Gaussian distribution? Why are the observation features (namely MFCCs) of an HMM-based recognition/synthesis system modeled by a Gaussian distribution?

Even the state duration is modeled by a Gaussian in this paper:

K. Tokuda et al., "Speech synthesis based on hidden Markov models", Proceedings of the IEEE, May 2013

For observations people usually use a mixture of gaussians, not a simple gaussian. They have few advantages - EM algorithm is fast to converge and GMM approximates wide variety of distributions pretty well. Probability of the model can be computed efficiently with gaussian selection. Last, GMMs are easy to cluster for context-dependent tree for phonetic context, which is important for LVCSR.

There are disadvantages as well, for example, GMM estimation is not very robust to the noise in features, also it does not account for some hidden parameters, for example "speaker" identity is not properly modeled with GMM.

There are extensions developed here, for example, subspace GMM is proposed to model speaker variability. In this model GMM parameters are computed as a projection of the GMM supervector thus allowing to model linear combinations of speakers and each speaker has individual GMMs.

Overall focus in state-of-art implementation has been shifted to DNNs for classification and they are better classifiers of course. Context-dependent tree to classify phonetic context is still constructed with GMMs though.


Modifying sizes and thickness of letters to "embed" an image in text I am wondering if this is defined, names and possibly solved - suppose I have a grayscale image and about one page worth of text. I would like to attempt to give  the text an appearance of an image when looked from far away. This would mean that text would be mapped to the space of an image and each letter would have size and thickness proportional to the mean pixel value of the area it maps onto. Any output  is ok (vector or raster).

Is this something that is used, named, does anyone know a place to look for this before I start writing the code from scratch? Python, R, MATLAB or a command line tool would all be acceptable.

If the whole pipeline is not a known use case, I would appreciate a relatively easy way to put letter X with (locX,locY,size,thickness) = (x,y,s,t) in an image, and I would write the rest myself.

(This question is too general for stackoverflow, but if it should be moved elsewhere I would appreciate that help as well.)

...give the text an appearance of an image when looked from far away.

The technique you are looking for is called "Halftoning".

It works with all sorts of tiles as long as you can get an indication of their "grayscale value" or, in the case of coloured images, the "grayscale value" for each colour (RGB).

Consider for example these two letters in a very elementary and rather retro "dot matrix" representation. Grey values represent background, White values represent foreground.

The amount of "dots" that each letter turns on is 10 and 5 respectively. Therefore, "A" would be brighter than "I" if they were to be used as pixels.

If you count how many "dots" each character (in a particular alphabet of characters you may have available) turns on and call that "brightness", you can then order the letters from darkest to brightest and create a correspondence between pixel value and letter that has to be typed.

Similarly, you can do this for any coloured image and create a colour palette (a set of RGB values, potentially ordered by their brightness).

This is the simple version of half-toning.

There is a tonne of software that do this (or slight variants) already, such as caca and hiptext  (or even img2txt).

The next step is to apply dithering which is a way of reducing the impact of the quantisation of the image from fine grain pixels to the tiles or maintain a list of alternative tiles for the same "grey value" to avoid repetition. Some of the tools above already cover these.

Beyond this simple method there are variants that given a set of tiles they will also try to match the "shape" of the tile to the shape of the image patch they are placed over. The output of those programs is really impressive. Unfortunately, I could not find the one I had in mind but there seems to be this one and one of the papers describing the technique is this one

Hope this helps.


Spectrogram power/magnitude normalization in analogy with image intensity normalization In analogy with the 8-bit grayscale image color depth, in which the intensity of a sample pixel is normalized to a value between 0 and 255 by taking into account the minimum and maximum intensities present in the original image, is there a standard for normalizing the powers/magnitudes (of short time Fourier transform buckets) in spectrograms of audio samples?

(I am new to audio signal processing, so please let me know if my question does not make sense.)

there isn’t any one solution to this problem.  the application governs the solution.

There are a few issues associated with mapping a DFT bin from a STFT calculation to a pixel on a display.

In many applications, the spectra tends to fall off with increasing frequency.  The low frequencies tend to dominate in the signal. Spectral leakage can cause these low frequencies to mask the high frequencies. Picking an appropriate window is important. Whitening your signal prior to taking DFTs is usually a good idea and can help picking a proper window. The window influences your overlap fraction.

One can also do some short term averages (of power) in time of sequential spectra to help improve SNR for some classes of signals. This also compresses spectra in time which saves display bits.

One issue is gray scale or colors? This influences how many pixel bits you have.  Color gives you more bits to map to but one needs to tune a color map.

Color is useful for overlays on gray scale like cursors.

Using the log of intensity (dB) does a very good job of mapping the dynamic range out a FFT to the range of a display pixel.

A fixed point FFT can be easier to map to pixels over a floating point FFT.

One reason that there isn’t a standard display technique is that display technology keeps advancing.  I first did this sort of display with a total of 480 bits across.

One trick is to set all the pixels below a threshold (like a median) to zero. This cuts down a lot of clutter. Scene complexity can cut either way.

I would recommend that making your spectra flatter over frequency is worth the effort.

One should also think about the interactions a user will have with the STFT image. Things like zoom, drill down, harmonic cursors, measurements like time differences and frequencies, ......


Estimate the Transfer Function of an Unknown System Suppose you have a system, H, that you want to estimate its transfer function.

You have a finite number of complex input samples, x, and noisy complex (magnitude and phase) output samples, y:

In practice, how do you estimate it's transfer function given a mean-squared error criteria?
What algorithms are used to estimate transfer functions when you have limited and noisy data?
How would you estimate the transfer function using a Wiener filter?

Here's the way I think about a discrete Wiener Filter

Consider a sequence of observations $\mathbf{y} \in \Re^n $

Form a matrix from the input $\mathbf{x} \in \Re^{n+r-1}$ by shifting columns one sample each:
$$
X= \begin{bmatrix}
x_1 & x_2 & ... & x_r \\
x_2 & x_3 & & x_{r+1} \\
x_3 & x_4 & & x_{r+2} \\
... & & & ...\\
x_n & x_{n+1} & ... & x_{r+n-1}\\
\end{bmatrix}
$$

Then the act of correlating the signal $\mathbf{x}$ with a filter $\mathbf{h}$ can be represented as $\mathbf{y} = \mathbf{X h}$

So $\mathbf{X^Ty}=\mathbf{X^TX~h}$ and the least squares solution for $\mathbf{h}=(\mathbf{X^TX})^{-1}\mathbf{X^Ty}$

Here's an octave/Matlab example that might help illustrate

r = 31; % filter length
n = 1e4;% signal length
SNR = 0;% AWGN measurement SNR in dB
x = randn(n,1); % random signal
h0 = randn(r,1); % unknown filter
y0 = filter(h0,1,x); % convolve x0,h0
y0 = circshift( y0, -(r-1)/2 ); % align x,y0
y = y0 + randn(size(y0)) * norm(y0)/sqrt(n) * 10^(-SNR/20); % add noise

% Wiener-Hopf solution
XX = xcorr(x,x,r-1);
Xy = xcorr(x,y,(r-1)/2 );
hrev = toeplitz(XX(r:end))  \ Xy; % least squares solution
h1 = conj(flipud(hrev));  % change from a correlating filter to a convolution filter

nmse_recovery = 20*log10( norm(h1-h0) / norm(h0));
plot([h0 h1]);
legend('truth','recovered')


How to extract the people and objects from an image separately? For example, I have an image

There are 7 people in the image and one apple and a bird.
How can I extract every person and every object from this image and put the 7 people and 2 objects into 9 different new images? How can I also identify which is person and which is object?

I think what you're looking for is called "category level object recognition". As far as I know, it's still an active research subject. Meaning you can't expect to find high-quality libraries that you can use out of the box without understanding the magic behind. And even state-of-the art research code only achieves around 90% correct classification.

A common approach is the so called "bag of (visual) words" model. Meaning: You take a lot of labeled training images containing various objects. You use a keypoint detector&descriptor like SIFT or SURF to find the "visual words", i.e. similar-looking patches in the training examples, and you train a machine learning algorithm (like a support vector machine) to recognize the category labels based on the presence or absence of these words in each training image.

Then you again use the same keypoint descriptor, e.g. using a sliding window-technique on your test image, find the closest visual words, feed that information to your learning algorithm and you get an object classification for the sliding window (if you're lucky, that is).

Sadly, if you want to do that with Mathematica, you'll have to reinvent the wheel a lot. The SURF keypoint detector/descriptor is included in the ImageKeypoints function, but I don't think you can make it use sliding windows. And there is no machine learning functionality to speak of.


Detect Gradual Increase Before a Decrease in Noisy Data I have data that'ss not too noisy and I am trying to detect a pattern where it gradually increases then decreases in a short period of time (20 ticks? it should be roughly similar per session but can vary slightly between different sessions). This pattern usually happens in pairs but not always.

So looking at the top plot of the data over time, the pattern I am after is at time 100 and 300. you can see the data is roughly stable, then gradually increases and decreases and stabilizes, then after a while again increases, decreases and stabilizes.
As you can see there are also instances where it suddenly changes (550 & 800), or it increases then decreases over a longer period of time (700-800) but I need to detect my short pattern only.
What I have done sofar:
I have taken the variance of the last 10 values which is the middle plot. It gives me a relatively good idea when the data starts to deviate but, variance ignores the sign so I cannot tell if the data is increasing or decreasing and it doesn't differentiate between sudden and gradual changes.
So I decided to take the squared difference between each value and the previous value and keep the sign. I took the square to exaggerate the difference to make it more distinguishable from smaller variations and kept the sign to help identify when it increases then decreases.
I can't say
if diff[i] > diff[i-1] for 5 ticks, and then diff[i] < diff[i-1] for 5 ticks because diff is noisy & changes its sign to positive midway at 130 even though overall it was decreasing.
I feel like there would be a better way than trying to filter the noise and then applying that if statement above. But I'm drawing blanks at the moment.
Whats the best way to programatically detect a short but gradual increase then decrease in noisy numbers?
I've added a csv with the data at this pastebin: csv_pastebin

What you need to do is find a feature which separate the phenomenon you're after from the noise.
If the jump is high relative to the noise STD it means it has good SNR when you compare the instantaneous derivative to the noise.
As you figured it out, in the case of great SNR as you shown above it works pretty easily and your approach did the work.
What about cases with bad SNR?
Then you need to apply some pre processing to increase the SNR.
You can apply Low Pass Filter but it means it will blur the exact location of the "Jump" you're after.
A nice trick I uses once was borrowing from Image Processing a class of filters called Edge Preserving Filters.
For instant it is really easy to apply the Bilateral Filter in 1D.
While it will remove noise and improve the SNR on flat regions it won't hurt the "Jumps" you're after.
So after using it even signals with worse SNR then the example above will start look like the case above and make your current method viable.


How to find the perfect regularization number? I have an issue!

Assume that we are going to solve $Ax=b$ but $x$ contains noise.To minimize that noise we can use regularization:

$$x = (A^TA + \lambda I)^{-1}A^Tb$$

Where $\lambda$ is a small number $> 0$, or sometimes $\lambda$ can be huge!

Here are three pictures of one situation where I solve $Ax=b$ with noise free and one $Ax=b$ with noise and without regularization and one with noise + regularization.

Noise free case: Here we can see that it works smooth. In this case, $x$ is a matrix because it contains two signals.

Here we are using 5% noise, that will cause much trouble. No regularization.

Here we are using 5% noise and a large number for regularization.

Question:

Is there a way to determine an optimized regularization parameter $\lambda$?

Source:
I made a library named Mataveid and it's made for system identification.
I have no problem with the noise actually, but I want to find a way to automize the choice of regularization parameter $\lambda$.

Approaches can be any of the following:

Model the noise that you expect, for ex: gaussian (least squares in Minimum variance unbiased estimator for a linear signal model in presence of gaussian noise). Based on this model try and estimate the noise variance, the regularization term should be close to noise varainace.
Deploy machine learning techniques based on signal parameters such as dynamic range, frqeuency, FFT noise floor etc. to estimate noise
Follow a recursive approach with adjusting the regularization to minimize a suitable objective, ex: gradient descent etc.

In general your regularization is as good as how close it is to the actual noise or estimation of noise. This goes back to the idea of whether a frequentist approach is better or Bayesian approach. The answers is if the modelling of the parameter is accurate to it's actual behaviour then the Bayesian approach is better but if you model a parameter with a distribution but it is actually not following that closely then Bayesian would lead to suboptimal result and better go for a frequentist approach.


FFT sample ordering with odd N I have series of $ N =133 $ data sampled symmetrically around the origin in an actual space at $ x(n \tau) = x_n $ with $ n \in \{ -(N-1)/2, \, ..., \, 0, \, ... , \, (N-1)/2 \} $. Now, I'm a bit unsure of how a typical FFT implementation performed on my data will be correct as a result of $ N $ being odd. My data will be zero padded to a length of $ (2N-1)_2 = 512 $, which is even. Ordering the sample from $ 0 $ to $ N $ as is done in an implementation,
$$
X_f = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-2 \pi i f n / N}
$$
I know that $ x_0 = 0 $, and $ X_0 \stackrel{!}{=} 0 $ needs to be my DC term, as $ f $ is sampled at intervals of $ m / N \tau, \, m \in \{0, \, ..., \, N-1 \}$. But now my data has one more element on one side of my zero-datum than the other since all the spatial frequencies above $ N/2 = 256 $ are negative frequencies due to the periodicity $ X_{-m} = X_{N-m} $: there are $ 256 $ points to the right of zero, and $ 255 $ points to the left (as $ X_{257} \rightarrow -255 $ cycles per length, ..., $ X_{511} \rightarrow -1 $ cycles per length).
So my real data is sample-symmetric about the origin, but a zero-padded data set isn't. It's hard to explain exactly what I'm having difficulty with, but it seems like my transform is going to have some kind of "fencepost" error in it due to some kind of offset by one in my data about the center. There will be one extra zero-pad data point on the right side than the left. In fact, if I had taken an even number of real samples from the start, there is no datum actually associated with the origin at all, so shouldn't an $ n = 0 $ term zero out some poor piece of data that doesn't deserve it? I know my fear is unjustified, but I can't seem to work through whatever imaginary problem it is I'm having.

Since your input is discrete the spectrum is periodic with the DFT length. $X_0$ is your DC term, but so is $X_{512}$, $X_{1024}$, $X_{-512}$, etc.
It really doesn't matter which exact period you pick: [0 511], [-255 256], [-256 255], etc. Any choice will work and you always get the same amount of data. Note that $X_{-256} = X_{256}$.


Help with algorithm for modulating oscillator pitch using LFO I am developing a software emulation of an analog synthesizer.  I am trying to modulate the pitch of an oscillator using an LFO.  For each sample fed to the computer's sound system, I am calculating the frequency to be input into the main oscillator like this (pseudocode):

osc_frequency = note_frequency * (1 + tuning) * (1 + lfo_y * lfo_mod_depth)

The variables in this statement are described as follows:

note_frequency = frequency of note to be played in Hz
tuning = oscillator fine-tuning in percent of played pitch (ex: -0.02 = detune by 2%)
lfo_y = current y value of lfo waveform (ranges from -1 to 1)
lfo_mod_depth = depth/intensity of effect to apply to oscillator in percent

However, this computation does not yield the desired result.  I expect to hear the pitch modulate up and down, locked around the center frequency (note being played).  What I am getting is a modulation effect that causes the pitch to "run away";  I can't tell exactly what is happening, but it sounds like one of these:

The modulation intensity increases over time (the high/low frequency mark reached by the modulation grows higher/lower the longer the note is held)
While the modulation intensity remains constant over time, the center frequency increases while the modulation oscillates around it

Am I using the correct approach or not?  If not, what should the correct approach be?  Any help with this is much appreciated.

As said by Jason, this might simply be that you are not implementing your oscillator correctly - for example by multiplying the frequency by time instead of integrating it.

Note also - and this is unrelated to your topic but really worth observing - that your formula for frequency modulation implements a behavior very different to that of most synthesizers, and will sound uncanny to a musician.

For example, if lfo_y oscillates between -1 and 1 ; lfo_mod_depth is 0.5 ; and if note_frequency is equal to 220 Hz, osc_frequency will sweep between 110 and 330 - that is to say between -1 octave and +1 fifth around the note. So the pitch modulation will appear to be centered in the hertz scale, but will not be centered in the perceptual musical scale.

The correct behaviour is to have something like:

osc_frequency = note_frequency * 2 ** (tuning / 1200.0 + lfo_y * lfo_mod_depth)

Then:

tuning is expressed in cents, a musically relevant unit (100
cents = 1 semitone).
your pitch modulation is "musically centered",
and lfo_mod_depth is expressed in octaves.

This is called "exponential FM", and this is the norm on synthesizers. On analog synthesizers, this is implemented by summing the LFO signal to the CV that hits the exponential converter of the VCO. On digital synthesizers, this is implemented by applying the modulations on a high-resolution internal representation of pitch which is still on a musical scale - before converting to a frequency or a phase increment.


Finding Laplace Transform without ROC While studying Laplace Transform i found that region of convergence (ROC) is important because for some problems we have same Laplace Transform but different ROC helps us to take correct inverse Laplace Transform.
So, Now i am practicing inverse Laplace Transform problem i found almost every probelm to find $x(t)$ is without ROC.
So i want to ask that without have ROC how can we solve inverse Laplace Transform?

Strictly speaking you can't because without specifying the ROC, the inverse Laplace transform is generally not unique. However, in many contexts there is the implicit assumption of causality of the corresponding time function (i.e., $x(t)=0$ for $t<0$), which is equivalent to stating that the ROC is a right half-plane.


Generating blue/violet noise without FFT I need noise with a high frequency bias for a realtime audio application, so I'm interested in things like blue noise or violet noise.

Is there a way to generate such noise without using a FFT or an IIR filter? In other words, something a little less computationally expensive?

Less expensive than filtered noise will be hard.

Practically all numerical implementations of colored noise that I've encountered so far are based on some random number generator designed to be as white as possible, followed by some noise shaping.

Violet noise itself is sometimes called differentiated white noise (according to wikipedia, at least), but implementing a differentiation in DSP is no different than designing a FIR with a specific high-pass characteristic, so that doesn't help you either.

In fact, that's what my "let's whip this up in less than a minute" attempt looks like; I'm a GNU Radio nerd, so that was my DSP framework of choice to get a blue noise generator done; I made the mistake of calling the files "violet"/"purple", although the noise is blue:

Left: visualization of signal spectrum (frequency scale was arbitrarily set to $\pi$ MHz); Right: GNU Radio Companion showing generating Flow graph.

To explain this flow graph: The idea was to shape the noise with a filter. Now, (FIR) filtering in time domain means convolution and that's terrible computationally. But convolution in time is multiplication in frequency domain. So, signal->FFT->element-wise multiply with constant vector->IFFT is usually faster than direct filtering, at least for filters of length >> 16 on PCs, as far as my experience goes.

Usually, one would use the existing FFT FIR block to implement that, which will take care of applying overlap/add, so that the filter response is actually identical to the time-domain FIR being replicated, but: we don't need that. Not doing the zero-padding before DFT'ing, and not doing the overlapping will not really change the white characteristics of the input signal.

Since the FFT blocks already come with the option to apply a window, I'm choosing a rectangular (i.e. constant 1) window for the FFT (a no-op, practically), and a upramp/downramp window (Python: range(0,fftlen/2)+range(fftlen/2,0,-1)) for the IFFT.

Note that albeit you said you wanted to use this for an audio application, I've done this with complex noise – the reason is simple: The (real-mode modified) FFT is less than twice as efficient for real data than it is for complex, so I just told myself "well, then you get two streams of purple noise; that's fine, too".

For your listening pleasure, I replaced the sink with a file sink, and converted the result to a stereo 24bit FLAC file; you can download the audio file, and find the GNU Radio companion file in that state in this Github Gist.

  a little less computationally expensive

Benchmarked this on my PC workstation, removing the graphical sink. Average throughput is ca. 10315 IFFTs per second, which, at an FFT length of 1024 that I used, equals to 10.56 MS/s (and, also, a lot of time is spent generating the white noise, because I just used the much-higher-quality-than-needed MT19937-based noise source with a Gaussian shape). That means that per second of processing time, you get two streams of purple noise worth 220s audio at 48 kHz. In other words, if you had e.g. a sound card throttling the amount of samples consumed, this flow graph, which effectively occupies two CPU cores, would be running only 1/220 of the time, ie. consume 2% of one core, if you were to restrict it to running on a single core.

I claim that to be computationally inexpensive.

Now, admittedly, if you're doing this in hardware, things might be different, because FFTs can be quite space-consuming – but it's probably still not what breaks your back at audio rates; by multiplying in frequency domain with a ramp, you're eliminating the need to store FIR coefficients (you'd just count samples and generate the factors ad hoc), and that's much cheaper than having a direct FIR implementation.


Filter for noisy accelerometer data I tried to read up and get a filter for a 3-axis accelerometer but this is more confusing than I thought.
I have an IMU with a 3-axis accelerometer on the end of a robot link which is rotated  around on one axis.
When at rest with just gravity acting on it I get these readings:

Now the band of the data is kinda broad and very noisy . I thought that at first I could maybe look at the frequency spectrum and cut out unwanted frequencies but honestly the FFT doesn't tell me anything I could remove. What does initial spike in the spectrum mean and and how would I put a high or low pass filter over the data if it looks like such a mess? Do I need to smooth it in some way?
Preferable I would want a nice arc of the data which follows the movement of the end of the arm but I don't know how to achieve that. Could it be that the arm movement by the motor introduces too much noise?
Thanks for any help!

Most accelerometers suffer from bias --- a non-zero, constant rest measurement. See, for example, this page.

  A further imperfection is called zero g bias level, or offset. This is the reading the sensor provides with no acceleration, and, unfortunately, that is rarely zero (in the Z direction it should be 1 g). In fact, the error can be substantial, frequently hundreds of times larger than the device’s resolution.

So you need to remove this DC offset (which is how the bias appears) from all your measurements.  If you need to know the $g$ value from the $z$ accelerometer, you will need to calibrate it first in the $x$ or $y$ plane (to remove the effect of gravity).

Have a read through this answer to another question, and try applying a DC blocker to your data.


Which image for DoG filter? I am studying computer vision and have some trouble with difference of Gaussian filter.

We apply Gaussian filters with different values of sigma on original image then subtracting from them.

However, Gaussian filter blurs image, so some edge might not be identified.

Why could we use DoG?
How effective is DoG?
Which type of images work best for this filter?

Why could we use DoG? An old and low-cost edge-enhancement technique consists in subtracting a blurred image from the original image. Its filter interpretation  is an impulse filter (the neutral filter for the original image) minus the blur filter. The resulting filter is positive at its center, and negative around, a sort of coarse Laplacian filter.

The above impulse filter can be seen as the limit of a Gaussian filter whose $\sigma$ tends to $0$.

The DoG filter somehow generalizes the described old edge-enhancement technique, with a pair of blurring filters. Perhaps surprisingly, two blurs, when subtracted, can sharpen edges. The resulting filter acts as a pass-band with appropriate $\sigma_k$ kernel ratios.

How effective is DoG? its effectiveness is close to that of the Laplacian of Gaussians, which it mimics, see below in 1D.

Which type of images work best for this filter? It is believed to work well for large images, since it has very fast implementations. Due to the "double filtering effect" in high frequencies, it is somewhat robust to noise. Depending on the kernel size choice, it enhances features that possess a certain scale extend.


What are typical audio waveform synchronization methods? I have two roughly similar (not identical) audio files, both being recordings of the same audio source material, but made independently of each other.
I need to calculate the time offset between each file, for instance one recording would be started at time A and the other at time B. By inspecting and processing the audio wave forms in each of the files I need to calculate the offset (A - B).
What are typical and effective ways of achieving this?
It can be assumed that both files are mono and that the sampling rate is consistent between the two files.

Crosscorrelation, find the maximum, and use that to offset one file vs the other.


Does/can physical audio theory exist? I've been wondering,

Since some sound generating physical phenomena have been physically modelled as physical modelling synthesisers,

Does there exist a "physical audio theory" (or scientific subdomains fit for belonging to "physical audio theory"), which is like a subfield of applied mathematics or physics that describes the theory of sound generators, propagation and such. For example, the Karplus-Strong algorithm and its analysis would belong there and it would be presented as a derivation of string physics. Sure there's the physics of sound under physics, but I have not seen a coherent mathematical theory on "physical modelling synthesis", even if "physical modelling synthesis" could perhaps be laid out as mathematical formalisms, e.g. as an "applied linear algebra" body of theory. It seems like much of audio DSP is just scattered information.

For those interested, I also think there are at least some books concerning "physical audio theory", an example would be this book:

https://mitpress.mit.edu/books/designing-sound

It does contain chapters that predominantly analyze natural phenomena or "sound generating settings" and try to translate those into expressions in mathematics, using concepts from physics. The derivations don't differ that much from e.g. typical engineering analysis. Much like the Karplus-Strong algorithm was developed by analyzing string instruments.

Physics-based sound synthesis is a pretty multidisciplinary topic.  Quoting from the preface of Physical Audio Signal Processing:

  [T]he material of this book is multidisciplinary, building on results
  from physics, musical acoustics, psychoacoustics, signal processing,
  control engineering, computer music, and computer science. Such
  diversity is typical of applied research.

It wouldn't really make sense to have it be some kind of subfield of math or physics; computer graphics is in a similar position.  Physics-based sound synthesis is mainly numerically simulating a physical model.  The trick, as usual, is making approximations to arrive at something reasonably efficient while not losing too much of what is important.  This is how fields like psychoacoustics and computer science come in.  There isn't a mathematical theory that will tell you what approximations won't produce unpleasant sounds and will run efficiently on a SHARC DSP.  The constraints are too ad-hoc.  Of course, you can put the constraints of a particular problem into a general mathematical theory of approximation and optimization to help derive solutions.


Faster way of getting 2D frequency amplitudes than DFT? I'm making a small program which gets the DFT of an image to get a general idea of the image's overall orientation.  It does this by rotating a line in a radar sweep type pattern, keeping track of which angle has the highest sum of frequency amplitudes. (It actually only does $0-180^\circ$ though because the testing line is symmetrical across the origin)

Anyways, I was wondering, is there a less computationally expensive way to get frequency amplitudes of a 2D image than doing a 2D DFT?

I know that it is separable so can do better than the naive $\mathcal O(N^4)$ operation - I think it becomes $\mathcal O(N^3)$.  I also know that the FFT can do a 1D DFT in $\mathcal O(N \log N)$ instead of $\mathcal O(N^2)$.  Lastly, I know that it's massively parallelizable so could be (massively) multithreaded on either CPU or GPU.

Those things aside, are there any less computationally expensive ways to get frequency amplitudes of a 2D image than to do a full 2D DFT?

Note that the computational complexity of a 2D FFT would be $\mathcal O(N^2\log N)$ instead of $\mathcal O(N^4)$ for the naive 2D DFT.

That said, since you mention that you want to obtain spectrum information along beams of a radar-like pattern, it sounds like you are really more interested in computing the 1D FFT of each such beam. In that case, for $M$ beams the computational complexity would be $\mathcal O(MN \log N)$.

As an illustration, lets process the following figure:

from which we shall extract a set of $M$ beams with the following sample matlab code:

% Number of beams to extract
M = 200;

% Convert the input image "img" to polar coordinates
c = size(img)/2+1;
N = max(size(img));
angles = 2*pi*[0:M-1]/M;
radius = [0:floor(N/2)-1];
imgpolar = zeros(length(radius), length(angles));
for ii=1:length(radius)
  xi = min(max(1, floor(c(2)+radius(ii)*cos(angles))), size(img,2));
  yi = min(max(1, floor(c(1)-radius(ii)*sin(angles))), size(img,1));
  imgpolar(ii,:) = img(sub2ind(size(img), yi, xi));
end
% Compute the FFT for each beam angle
ImgFD = fft(imgpolar,[],1);

figure(1);
freqs = [0:size(ImgFD,1)-1]/size(ImgFD,1);
surf(angles, freqs, 10*log10(abs(ImgFD)+1), 'EdgeColor', 'None');
view(2);
colormap("gray");
xlabel('Beam angle (radians)');
ylabel('Normalized frequency');

to yield:

which can be collapsed to the sum of amplitudes as a function of the beam angles to give:

SumAmplitudes = sum(abs(ImgFD),1);
figure(2);
hold off; plot(angles, 10*log10(SumAmplitudes+1));
xlabel('beam angle (radians)');
ylabel('Sum of amplitudes (dB)');

As a side note, if you can use sum of squared amplitudes along those beam (instead of sum of amplitudes), then you can do it directly in spatial domain thanks to Parseval's theorem (which would bring the computational complexity down to $\mathcal O(MN)$, dominated by the conversion to polar coordinates). The equivalence (for the sum of squared amplitudes) can be seen using the following code:

% Compare the result of square amplitude summation in the frequency domain vs spatial domain
SumFD = sum(abs(ImgFD).^2,1)/size(ImgFD,1);
SumSD = sum(abs(imgpolar).^2,1);

figure(3);
hold off; plot(angles, 10*log10(SumFD+1), 'b');
hold on;  plot(angles, 10*log10(SumSD+1), 'r:');
xlabel('beam angle (radians)');
ylabel('Sum of squared amplitudes (dB)');
legend('Frequency domain', 'Spatial domain', "location", "southwest");

Notice the overlap of the curves computed in the spatial and frequency domains:

Update:

If you are in fact computing beams of a radar-like pattern in a 2D frequency plot as you seem to suggest in this other post, then you best bet comes back to doing a 2D FFT which would be order $\mathcal O(N^2 \log N)$. You could then perform the conversion to polar form for the sum of frequency coefficients which would add a small $\mathcal O(N^2)$, so the result is still dominated by the 2D FFT.


Performance of MATLAB's `conv2()` vs. `imfilter()` A lot of people use imfilter to achieve a 2-D convolution between an image and a filter, but the majority of people use conv2 instead of imfilter because it is faster than imfilter by at least 20%.

Examples:

filtered_image_with_imfilter = imfilter(image,filter,'symmetric','same','corr');
filtered_image_with_conv2 = conv2(image,filter,'same');

Could someone explain to me in detail why conv2 is faster than imfilter.

First, I think conv2 and imfilter will give you the same result if you change the filter option of imfilter to conv instead of corr. imfilter uses correlation to filter images by default that starts from one side of the image, whereas covolution starts from the other, so there may be some small differences in the filter output.

Secondly, my test shows imfilter on 2D image is faster if the filter size if the same as images.

clc
data = (rand(200));
kernel = ones(200)/4e4;
tic
blurredData = conv2(data, kernel,'same');
conv2t=toc
tic
blurredData = imfilter(data, kernel,'same','conv');
imfiltime=toc

>> conv2t =

       2.7425

  imfiltime =

       0.1552

Both imfilter and conv2 are hardware optimized and multi-threaded. For floating point inputs, the performance of them are supposed to be roughly the same. But in most cases in image processing, the filter kernel is linearly separable, imfilter does separation of the kernel into row and column vectors. As a result, they are much faster.  Yet if the filter kernel is (1) not linearly separable, or (2) a small kernel such as a 5*5 filter kernel. conv2 wins in the processing time. I guess the reason is some time is spent on checking the separability of the kernel in imfilter. I am not sure whether there are other different mechanisms that leads to the processing time difference between the two.

Finally, if the filter size is the same as image. fft is always the best one (O(nlogn)):

clc
data = (rand(200));
kernel = ones(200)/4e4;
tic
blurredData = conv2(data, kernel,'same');
conv2t=toc
tic
blurredData = imfilter(data, kernel,'same','conv');
imfiltime=toc
tic
D=fft2(data);
K=fft2(kernel);
blurdata=ifft2(D.*K);
fftime=toc

>>conv2t =

          2.7429

  imfiltime =

          0.1327

  fftime =

          0.0164

EDIT

@geometrikal's detailed analysis:
The reason fft is faster than imfilter with same size images is that imfilter pads the inputs (even with 'same'). For small filter kernels, not much padding is needed for imfilter, whereas to use the fft the kernel will have to be padded to the size of the image. In these cases imfilter is faster.


Multilevel partial wavelet reconstruction with PyWavelets I'm looking for a way to partially reconstruct branches of a wavelet decomposition, such that the sum would recreate the original signal. This could be achieved in MATLAB using:

DATA = [0,1,2,3,4,5,6,7,8,9]
N_LEVELS = 2;
WAVELET_NAME = 'db4';
[C,L] = wavedec(DATA, N_LEVELS, WAVELET_NAME);
A2 = wrcoef('a', C, L, WAVELET_NAME, 2);
D2 = wrcoef('d', C, L, WAVELET_NAME, 2);
D1 = wrcoef('d', C, L, WAVELET_NAME, 1);
A2+D2+D1

ans =

    0.0000    1.0000    2.0000    3.0000    4.0000    5.0000    6.0000    7.0000    8.0000    9.0000

I'd like to achieve the same using pywt, but I'm not sure how to go about this.
The pywt.waverec function creates a full reconstruction, but doesn't have a level parameter for partial reconstruction.
The pywt.upcoef function does what I need for a single level but I'm not sure how to expand this for multiple levels:

>>> import pywt
>>> data = [1,2,3,4,5,6]
>>> (cA, cD) = pywt.dwt(data, 'db2', 'smooth')
>>> n = len(data)
>>> pywt.upcoef('a', cA, 'db2', take=n) + pywt.upcoef('d', cD, 'db2', take=n)
array([ 1.,  2.,  3.,  4.,  5.,  6.])

I managed to write my own version of the wrcoef function which appears to work as expected:

import pywt
import pdb
import numpy as np

def wrcoef(X, coef_type, coeffs, wavename, level):
    N = np.array(X).size
    a, ds = coeffs[0], list(reversed(coeffs[1:]))

    if coef_type =='a':
        return pywt.upcoef('a', a, wavename, level=level)[:N]
    elif coef_type == 'd':
        return pywt.upcoef('d', ds[level-1], wavename, level=level)[:N]
    else:
        raise ValueError("Invalid coefficient type: {}".format(coef_type))

level = 4
X = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
coeffs = pywt.wavedec(X, 'db1', level=level)
A4 = wrcoef(X, 'a', coeffs, 'db1', level)
D4 = wrcoef(X, 'd', coeffs, 'db1', level)
D3 = wrcoef(X, 'd', coeffs, 'db1', 3)
D2 = wrcoef(X, 'd', coeffs, 'db1', 2)
D1 = wrcoef(X, 'd', coeffs, 'db1', 1)
print A4 + D4 + D3 + D2 + D1

# Results:
[  9.99200722e-16   1.00000000e+00   2.00000000e+00   3.00000000e+00
   4.00000000e+00   5.00000000e+00   6.00000000e+00   7.00000000e+00
   8.00000000e+00   9.00000000e+00   1.00000000e+01   1.10000000e+01
   1.20000000e+01   1.30000000e+01   1.40000000e+01   1.50000000e+01
   1.60000000e+01   1.70000000e+01]


Ultrasonic signal detection I have created a fairly simple TDOA system that uses ultrasonic signals emitted from two speakers to geolocate (relative to the speakers) mobile phones.  The two signals are separated by frequency.

The system has the following constraints:

The signals must be inaudible.  To that end we stick to frequencies above 17 kHz.  A few people can still hear that, but most can't.
Sample rate is 44.1 kHz.
Music will typically be playing, so there is lots of interference at the lower frequencies.
We don't have control over how well the speakers and microphones work at the upper frequencies, so we've kept our upper limit at around 20 kHz.

The particular signal that I am using is BPSK modulated 13-bit Barker codes because of their good autocorrelation properties.  The autocorrelation looks like the following-

When I cross-correlate the expected signal against the received signal in real life, though, what I get typically looks like this-

The blue is the cross correlation with the speaker 1 signal, and the red is the cross-correlation with the speaker 2 signal.  It appears that the echoes are significant and, unfortunately, often stronger than the direct path signal due to the directional gain of the microphone.

I tried simply detecting the earliest appearance of the signal as that is likely to be the direct path.  This approach is very sensitive to the threshold that I use for deciding when the signal is present and so is not robust at all.

I would like a robust approach for determining the "true" arrival time of the signal- i.e. the arrival time of the direct path signal.  Perhaps some form of channel estimation and deconvolution?  If so, how would that work?

Data/Code: I want to make it clear that I am not expecting anyone to analyze the data or inspect my code.  I have made them available in case you want to do so.  I am mostly interested in ideas.

I made the raw received signal and modulated expected signals available for download.  They are all sampled at 44.1 kHz.  Correlating the received signal with the expected signals will produce something similar but not identical to the picture above because I move the received signals to baseband and decimate before correlating with the expected signals.

Received signal

Expected signal #1

Expected signal #2

Matlab scripts
The Matlab scripts has both the signal generation script (genLocationSig.m) and my receive/processing script (calcTimingOffset.m).

These are not the codes you are looking for...

As I mentioned in the comments, there are quite a number of ways to do robust TDOA. (Cross-correlation with Linear Chirps, Exponential Chirps, and CDMA-type methods). You have already built a TDOA system utilizing codes, (and that is indeed a good choice over linear-chirps if you need robustness to doppler), however you are limiting yourself artificially in two ways:

Barker codes only go up to length $13$. We can however make PN-sequence codes of arbitrary length to get a lot more coding gain.
The use of only $1$ bit in your transmission. We can encode an entire preamble of many bits to transmit, gaining further resilience to multipath.

Use a PN-Sequence:

Thus, very simply, change the codes you use to modulate your carrier by: Use PN-Sequences instead. PN generated codes can be of (nearly) arbitrary length, and can be generated via LFSRs. (They also go by the name 'whiteners' in some texts). Here are three PN-sequences of length $31$, $61$, and $127$ respectively.

PN_31 = [ 1  1 -1 -1  1  1 -1  1 -1 -1  1 -1 -1 -1 -1  1 -1  1 -1  1  1  1 -1  1  1 -1 -1 -1  1  1  1];

PN_61 = [ 1  1  1 -1  1  1 -1  1 -1 -1  1 -1 -1  1  1  1 -1 -1 -1  1 -1  1  1  1  1 -1 -1  1 ...
     -1  1 -1 -1 -1  1  1 -1 -1 -1 -1  1 -1 -1 -1 -1 -1  1  1  1  1  1  1 -1  1 -1  1 -1 ...
      1  1 -1 -1  1  1 -1];

PN_127 = [-1     1     1     1    -1     1    -1    -1     1    -1     1     1    -1    -1    -1     1     1    -1     1     1     1     1    -1     1     1    -1     1    -1 ...
       1     1    -1     1     1    -1    -1     1    -1    -1     1    -1    -1    -1     1     1     1    -1    -1    -1    -1     1    -1     1     1     1     1     1 ...
      -1    -1     1    -1     1    -1     1     1     1    -1    -1     1     1    -1     1    -1    -1    -1     1    -1    -1     1     1     1     1    -1    -1    -1 ...
       1    -1     1    -1    -1    -1    -1     1     1    -1    -1    -1    -1    -1     1    -1    -1    -1    -1    -1    -1     1     1     1     1     1     1     1 ...
      -1     1    -1     1    -1     1    -1    -1     1     1    -1    -1     1     1     1];

The circular and linear auto-correlations of the sequences are shown below. They will clearly yield white spectra, but more than that, we are no longer limited to $13$ chip lengths. In fact, the last code, PN_127, yields a coding gain of $10 \ log [\frac{127}{13} ] \approx 10$ dB gain over the barker sequence, all the while guaranteeing white spectra.

Transmit a preamble:

In your particular application, you mentioned that you were only transmitting one bit. You should try to avoid this if you can help it, and transmit as many bits as your application can allow, to get further coding gain out.

This is what is commonly done on communication protocols to align with the beginning of a packet. A (known) preamble is transmitted, composed of many bits. Each bit, is composed of many chips. (In our example, $31$, $61$, or $127$ chips with either of the above PN codes). Lastly, the bit sequence itself can be composed of yet another PN sequence, or if you like, you may transmit $13$ bits composing a barker pattern, with each bit being composed of either one of the the above PN sequences.

Try one or both of those solutions, and put up your results. I expect there to be tangible improvements that we can then iterate on. (Pulse shaping, different/longer PN sequences, etc).


Filter order and cut off frequency How can I know exactly what frequencies are not influenced by the filter?

I mean, I have a signal in the range of frequencies $5-100 Hz$ and I apply a high-pass, $4^{th}$ order Butterworth filter with cut-off frequency of $8 Hz$. The filter will not cut sharply at $8 Hz$, instead it will have a slope related to the filter order: the higher the order, the steeper the slope right?.

So, how can I compute the frequency from which the signal is left unchanged (so no attenuation due to the filter)?

Is there a way to obatin this info directly in matlab?

This answer is written in terms of low-pass but  the ideas directly translate to high pass.

A butterworth filter (low pass) of order $N$ with cutoff $\Omega_c$ is one which has a maximally flat passband, and has the squared magnitude response $|H(j \Omega)|^2 = \frac{1}{1+\left(\frac{j \Omega}{j \Omega_{cutoff}}\right)^{2N}}$. [I'm just writing it this way so that if you try to get a realizable transfer function via spectral factorization it is easy].

So, you specify the cutoff $\Omega_c$ and it gets sharper as you increase the filter order. Strictly speaking, from the the squared magnitude response, you see that $|H(j \Omega)|^2=1$ if and only if $\Omega=0$, i.e. only the DC component is un-modified. However, from that formula, you can calculate things like frequencies where the magnitude is within say 99% of the original magnitude, by finding $\Omega$ such that $|H(j \Omega)|^2 \geq 0.99^2$ (note that $|H(j \Omega)|^2 $ is monotone decreasing in $\Omega$, so your range will be of the form $[-\Omega_{0.99},\Omega_{0.99}]$). Similar properties of other filters will hold as well - see Appendix B of Discrete-Time Signal Processing 2e by Oppenheim, Schafer and Buck for more details. Make sure you're familiar with the basic types of filters as well -- elliptical, butterworth and chebyshev to make sure your characteristics you desire are matched to what you want.

If you want to design filters with particular passband and stopband characteristics in the FIR case (such as controlling the maximum error), you may want to look at the Parks-Mclellan algorithm (described in Discrete-Time Signal Processing 2e by Oppenheim Schafer and Buck, ch. 7).

As for Matlab usage, you can plot the magnitude response of the filter and determine these things. It can be done from the filter visualization tool. TI also makes a filter design tool which is quite nice, called FilterPro.


When did signal processing become a field of its own? This is perhaps a bit soft question but when did signal processing become a field of its own? It would appear to have first been a branch of mathematics and statistics but then emerged as an independent field.

If there is any book that would cover this, or perhaps some special issue article, I would love to read it.

PS. Please help to add relevant tags; I could not find any

There's a nice historical account given by the IEEE Signal Processing Society, that references 1948 as the year the field began. I'm not sure if that's necessarily accepted as "official" by other sources, but their reasoning is pretty convincing (at least as a ballpark):

  The year 1948 may be regarded as the annus mirabilis in the emergence
  of the discipline of signal processing. For in that year Claude
  Shannon published the epoch-making "A mathematical theory of
  communication"; Bernard Oliver, John Pierce, and Claude Shannon
  published the classic argument for the use of pulse-code modulation;
  modern digital methods of spectrum estimation were introduced;
  error-correcting codes were introduced; audio engineering achieved a
  new prominence; and the IEEE Signal Processing Society, albeit under a
  different name, was established.

That quote is from the beginning of the second chapter, Annus Mirabilis: 1948.


how to extract 45 cycles from bunch of waveforms I want to extract 45 number of cycles from the bunch of cycle which is similar to sine wave of varying amplitude. please help me to do so. in this picture i have capture only 45 cycles but i have continuous wave forms.i.e after each 45 cycles the signal is repeating.

thanking you,

For best performance in the presence of noise, I recommend "locking" onto the cyclical signal using a Phase-Lock-Loop (PLL) approach. Then you can feed the PLL output (which will be a cleaner representation of your signal and can be a digital square wave if desired) into a digital counter, which you can then count to 45 after resetting the counter at any position desired.


Convert a 1D array into a 'rotated' 2D matrix simply put, let's say I have this array - [1 2 1]
I want to try and get this output in Matlab

[1 1 1
1 2 1
1 1 1]

It's like my array was 'rotated' into forming the above matrix such that if I start at any peripheral element and move through the center of the matrix (in a straight line of course) towards the other end, then I get back my original array ( [1 2 1] )
Can anyone help me achieve this?

I assume that the length of the array is odd and that the array is symmetrical. Also note that from your requirements not all elements of the resulting matrix are specified. Those elements will be zero in the code shown below.

v=[1,2,3,2,1];      % odd-length array
A=diag(v);
A=A+rot90(A);
n=(length(v)+1)/2;
A(n,:)=v;
A(:,n)=v;
A =

1   0   1   0   1
0   2   2   2   0
1   2   3   2   1
0   2   2   2   0
1   0   1   0   1


Deconvolve a symmetric sequence into minimum-phase and maximum-phase sequences I learned in school that any real-valued, even-symmetric sequence can be expressed as the convolution of a minimum-phase sequence and its time-reverse, which is maximum-phase.  In my case, the real-valued, even-symmetric sequence is $sinc(x)$.  What causal, real-valued, minimum-phase sequence convolved with its time-reverse produces $sinc(x)$?

I am familiar with the traditional approach of computing the frequency-domain zeros and constructing the minimum-phase sequence from those inside the unit circle and half of the pairs on the unit circle.  However, in this case, the sequence has infinite length.  I would prefer a closed-form expression for the minimum-phase sequence function.

In practice, I will be using a windowed version of the sequence -- an FIR filter.  I would also like to learn how to deconvolve any real-valued, even-symmetric, finite sequence (FIR filter) into its minimum-phase component, whose frequency response magnitude is the square-root of the symmetric FIR.  I would like a time-domain method for doing this because there are precision problems with computing frequency domain zeros in very long sequences, such as a million taps.

For approximating a minimum phase sequence from a finite length linear phase sequence (such as a windowed Sinc), here are some existing stackoverflow and dsp.stackexchange answers regarding a matlab minimum phase approximation script which uses cepstrum deconvolution:

https://stackoverflow.com/questions/16753508/cepstrum-deconvolution-matlab-code

Derive minimum phase from magnitude

And here's an earlier post on the same matlab script, with more comments, from the music-dsp mailing list:

http://www.music.columbia.edu/pipermail/music-dsp/2004-february/059372.html

And the code:

x = [x; zeros(size(x))];
m = size(x,2);
n = size(x,1);
wn = [ones(1,m); 2*ones(n/2-1,m); ones(1,m);  zeros(n/2-1,m)];
y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x)))))))));
y = y(1:n/2;);

I've found that a greater amount of zero=padding of the original linear phase sequence may help improve the approximation.


Why I don't get the right PSD I need to model a noise with a given PSD. To do this, I am starting from a white gaussian noise (WGN) and feed with the WGN a transfer function, which will act like a filter. In fact,it's easy to prove that if you choose a PSD of the white noise equal to $1 {\rm unit}^2/{\rm Hz}$,then the output is $ H(f) \times 1$.
Thus
$$ PSD =|H(f)|^2 $$

Gain in the transfer function acts properly, translating the PSD up or down.
The problem is with the zeros and poles. No matter how I create the transfer function (1st, 2nd, or higher order), the PSD will be flat in the region of interest $f_{\rm range} = [10^{-6}, 10] {\rm Hz} $

The transfer function is
$$ H = \frac{10^{-11}(s+2 \times 10^{-3})^2}{(s+10^{-9})^2}; $$

It is easy to build in Simulink band limited white noise with PSD = 1 (choosing sample and power).

Evaluate then the output of the TF in matlab with this

%% PSD estimation
for i =[1e-8 1e-7 1e-6 1e-5, 1e-4, 1e-3, 1e-2 1e-1 1]
    f_range = linspace(i,i*10,1000);

    [pxx_white] = pwelch(white,[],[],f_range);

    [pxx_color] = pwelch(color,[],[],f_range);

    loglog(f_range,sqrt(pxx_white./(2*pi)),'k')
    hold on
    loglog(f_range,sqrt(pxx_color./(2*pi)),'r')

end

xlabel('Frequency (Hz)')
ylabel('RPSD (unit/Hz^(1/2))')
legend('white','color')
grid on

Now my question is: Am I doing something theoretically wrong? Maybe I can't reach PSD in the low frequency region due to the simulation time? I simulate for 10000 s (1e-4 Hz).

In black: white noise PSD, in red: flat the colored noise PSD, the curve is the plot of the tf.

Since you said you wanted the frequency range to be from 1 $\mu$Hz to 10 Hz, I ran my simulations with 20 Hz sampling frequency and did 100 32k FFT PSDs.  My results looked like what you have. So I lowered the sampling frequency to 0.2 Hz. The sims started at 0 s, ended at 163835 s, with $\Delta$t = 5 s. The WGN was simply unit normal. My results are shown below:

Hope this helps.


What does filtering do to a signal in time domain? Does the signal gets distorted? When considering a transceiver system, if we do a certain filtering operation on the signal, will the BER of the system be affected? My point is, since filtering in frequency domain filters a portion of the spectrum, whereas in time domain it is basically a series of MAC operations depending on the number of filter taps. So, eventually my original signal after filtering should have a totally different waveform than the original. Now, if I want to get the original signal back in an ideal channel condition, do we have to do some sort of inverse filtering operations that compensate the operations done in the time domain?

I will be grateful if you could help me in this regard. Thank you in advance!

It's always important to keep in mind that more than 80% of signal processing is filtering. However, not every filtering operation should be "viewed" in the same way.

Let us take the example you mention, of an "ideal known channel" (assuming it's not in deep fade). In wireless communications that would mean perfect knowledge of channel to the receiver. Typically, you could invert the channel in a very basic and often used scheme known as Zero forcing i.e. if $$y=hx + n$$ where y is received signal then the equalization is simply $y/h$ assuming real signals and noise. Similarly an MMSE equalizer ( optimal for most SNRs) has a typical mathematical expression that you would find in standard textbooks. The point is it would be less intuitive to look at this equalizer from the perspective of it's frequency response. It might be more intuitive to look at it from the perspective of vector spaces in this case. Especially as the degree of freedom of a network grow. For ex: a MIMO system.

Filtering is essential to iradicate BER, simply because if you don't filter ( as mentioned above), the time domain signal has been essentially "changed" by the channel.

We are always trying "mitigate" the effects of the channel. I would say "mitigate" rather than "invert" because it is sometimes not even possible or doesnt make sense to simply invert. For ex: in the scheme above $y/h$, notice that the noise is also devided by $h$. If $h$ is small, then you boost the noise as well. So in such a case a MMSE equalizer would make sense. When viewed in line of probability theory, the MMSE equalizer is calculating the expected value of the posterior distribution of $x$ (your signal) after observing $y$. So in that sense it is not inverting the channel.

Long story short, we try to do what's best for the system depending on bit rates, how close the bit rate is to channel capacity, SNR, modulation schemes etc. Rather than simply invert the channel.

There are other scenarios such as an antialsiaing filter, in such cases it is best viewed from it's frequency response, like cut off frqeuency and phase responses. You would get a better idea of the time domain filtered sequence looking at it from the frequency domain perspective. For ex: the cut off frqeuency would tell us what maximum gradients are tolerated in the input signal to the filter and hence are part of the filtered sequence or are smoothen out. The phase response would tell us which freqencies of the filtered sequence, would be delayed by what amount. Whether delay would be constant for the whole signal (linear phase) or there would be different group delays (non linear phase).

Indeed the filtering will change the shape of the input signal, otherwise usually nothing interesting is happening, but the filter itslef is not designed to "distort" the signal, it is designed to "better" the signal according to application, indeed the shape /magnitude of time domain filtered signal will change after filtering, but for the better of the system.


Digital low pass filter vs continuous low pass filter,? What is difference between both implementations of low pass filter? Discrete/digital implementation and continuous implementation?

As far i am able to understand

In case of discrete implementation, passband spectrum exists between "-pi "and "pi"

While in case of continous implementation, passband spectrum exists between 0 and "f" . Where "f" is upper cut off frequency

A discrete spectrum has a unique span extending from $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate (this is referred to as the first Nyquist Zone in the analog domain), thus in terms of normalized frequency when we divide the analog frequency by the sampling rate, it extends from $-0.5$ to $+0.5$ or $\pi$ to $\pi$ depending on which units we use as described further below.

In the case of a discrete implementation, the unique spectrum extends from $-\pi$ to $+\pi$ radians/sample. This is a common unit of frequency for sampled systems while cycles/sample is also used, which extends from $-.5$ to $+0.5$. The discrete frequency spectrum is periodic beyond these boundaries, so there is no reason to include them as no further information will be provided.

These units are the units of "normalized frequency" as they are arrived at by dividing the analog frequency units (either cycles/second or radians/second) by the sampling rate $f_s$. Cycles/Sec is equivalent to Hertz. The sampling rate is in units of samples/second therefore we get:

Normalized frequency:    Cycles/Sec / Samples/Sec = Cycles/ Sample.

Normalized Angular Frequency: Radians/Sec/ Samples/Sec = Radians/ Sample.

That said, the passband of a discrete low pass filter must extend to an amount less the $\pi$ (in radians/sample units) since $\pi$ represents the highest frequency of the discrete spectrum.


Harmonics in 2-D FFT Signal In case of 1-D FFT it is far more easier to find number of harmonics and their Respective frequencies and amplitudes.
Lets say if i have 8 points signal the harmonics of it then would be 7 with one DC point i.e,

the harmonics after 4th points would have negative imaginary part as it below the nyquist rate.

but if I've 2D signal of 8x8 Matrix there would be total of 64 harmonics and the Nyquist rate will then be on the 32nd point. Am i right ?

If my above statement is correct, the fractional frequency for the first harmonic then would be ,

0 -> Dc
1st Harmonic - > 1/64 -> 64 samples in 1 cycle

Number of Cycle 1.
Number of samples = 64
Sampling rate = 0.015625

What would then be the actual frequency of this Harmonic , I am asking this because i couldn't be able to identify it when i am plotting it on Matlab.

Your diagram is incorrect. Computing the 2D-FFT of a 8x8 matrix is not the same as computing the 1D-FFT of a 1x64 vector. A 2D-FFT is "probing" your signal with templates which have both an horizontal and vertical frequency. A signal with a low horizontal frequency and a null vertical frequency is a progressive horizontal gradient. A signal whose vertical and horizontal frequency are both the nyquist frequency is a checkerboard pattern...

If you take an image and perform a 2D-FFT on it, you will find in the resulting matrix the following elements:

Cell (1, 1): DC component (average of all values in the matrix).
Cell (2, 1): first horizontal harmonic; frequency 1/8.
Cell (3, 1): second horizontal harmonic; frequency 2/8.
Cell (4, 1): third horizontal harmonic; frequency 3/8.
Cell (5, 1): fourth horizontal harmonic; frequency 4/8 (Nyquist).
Cell (6, 1): third horizontal negative freq. harmonic; frequency -3/8.
Cell (7, 1): second horizontal negative freq. harmonic; frequency -2/8.
Cell (8, 1): first horizontal negative freq. harmonic; frequency -1/8.
Cell (1, 2): first vertical harmonic; frequency 1/8.
Cell (1, 3): second vertical harmonic; frequency 2/8.
...
Cell (5, 5): fourth horizontal/vertical harmonic (Checkerboard pattern).
...

A good way to visualize this is to compute the 2D IFFT of a single pixel at each of the 64 positions of the matrix, and look at the real part.