texts: Impulse response of a causal system from transfer function in z-domain

texts:
The transfer function is $$H(z)=(z+1)/(z^2+z+0.5)$$
I need to find the impulse response h[n] of a causal system with x[n] as unit impulse.

I have tried to find the impulse response by the following methods:

partial fraction decomposition and then inverse z transform of transfer function.
homogeneous solution of differential equation.
Matlab iztransf() function
I get different solutions from the above methods and i don't know how to remove iota from the impulse response.What is the correct procedure to find the impulse response?


texts: Calculating the fundamental frequency using Autocorrelation often gives half the expected value

texts:
I'm currently writing a mobile app which needs to analyse musical notes and find the fundamental frequency to determine the pitch. To do this I'm reading in audio data, taking an FFT, taking the auto-correlation of that FFT and then finding the IFFT of that. Having done that I then find the maximum peak (which seems to always be a zero) and then find the next highest peak and determine the lag between them. Finally I calculate the period from this and find the fundamental frequency.

My problem is that this method is almost always returning a value which is half the expected fundamental frequency, only occasionally returning the correct value. Here is a graph of an autocorrelation that returned the correct fundamental frequency value:

When zoomed in on the relevant part we can see that there's a peak at around 370 which gives for the correct fundamental frequency of 117hz.

Here's a graph which returned the half the fundamental frequency value:

Zooming in on the relevant part we can see there's a trough at around 370 where we'd expect the peak to be:

Does anyone know why this happens and how I can resolve the problem/any other methods I can use?


texts: Detecting Zero Crossing in Image Filtered by Laplacian Operator

texts:
This is my first question on this site, complete noob. But here goes.

I have a positive Laplacian operator

[[0,1,0],
 [1,-4,1],
 [0,1,0]]

Now this Laplacian operator is used to find the outward edges of an image , IIRC. The image which I have in binary format is this

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Its hard to make out here, but the image represents some staircase type of  image. I have a PY script which does the convolve and gives me back the following result of convolution :

[[ 0.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  1.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  1.  0.]
 [ 0.  0.  0.  1.  0.  0.  0.  0.  1.  0.]
 [ 0.  0.  1.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  1.  0.  0.  0.  0.  0.  0.  1.  0.]
 [ 1.  0.  0.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Drawing straight line through all 1s actually does give me the outward edge of the binary image.

However, the zero crossing algorithm has confused me .
According to this link , the algorithm to detect zero crossing is:

  Do the following for each pixel $I(u, v)$:

  Look at your four neighbors, left, right, up and down
  If they all have the same sign as you, then you are not a zero crossing
  Else, if you have the smallest absolute value compared to your neighbors with opposite sign, then you are a zero crossing

But here, the 0 and 1 have the same signs , right ? Both are positive. I am confused in this.
If someone can take a small sample out of the result image, and mark some zero crossing for me, it will be clear to me.

Thank you in advance, it's a long post. I hope its clear though.


texts: Circular wrapping of an asymmetric function (in DFT calculations)

texts:
Convoluting a signal (using discrete FT) for a given interval [a, b] with a  Gaussian can be done by circular wrapping as shown in Numerical Recipes. I found a shortcut for circular wrapping so that the phase remains zero for an even-N data. Consequently, the peaks in the signal will not shift their positions after convolution. The shortcut is to center the Gaussian at the mean (a,b)+dt, where dt is a sampling interval. After that, one can do a circular shift by N/2.
Consider a hypothetical case where we have a skewed Gaussian instead of a Gaussian. Should I have the centroid at the mean (a,b)+dt before circular shift by N/2?
Continuous version

Circularly wrapped version


texts:
I'm struggling to filter some data from an accelerometer. I need to do it in python (and not just on the device collecting the data) to fix the data I have, so that I can build some classification algorithm

I have the following data on ROLL from an accelerometer (1st graph). I tried the following:

transform the outliers using winsorize (2nd graph)
and then, filter the data using Butterworth filter (3rd graph)

Problem is that winsorization technique does not fit my problem here, it creates points at the top and the bottom for the extreme values, and it generates those peaks after filtering. And it does not take into account when those outliers are within the given limits as well.

would someone have a solution to correctly remove those outliers? Any python function?
I'm happy to remove completely those outliers, rather than transform them. Below is the dream, expected output after filtering:

If I focus on 1 piece of outliers, we can see the following (my data distribution is a bit weird, I have a couple seconds every few seconds):

  jupyter notebook below

%pylab inline
warnings.filterwarnings("ignore",category=DeprecationWarning)
import pandas as pd, numpy as np        # Data manipulation
import matplotlib.pyplot as plt
from scipy.stats.mstats import winsorize
import scipy.signal as signal

df=pd.DataFrame.from_csv("/tmp/file.csv")
df.reset_index(level=0, inplace=True)
df=df.iloc[2000:18000]

df_with_outliers_trimmed=df.copy()
df_with_outliers_trimmed["ROLL"]=winsorize(df_with_outliers_trimmed["ROLL"], limits=0.05)

df_filtered=df_with_outliers_trimmed.copy()
# design the Butterworth filter
N  = 1    # Filter order
Wn = 0.03 # Cutoff frequency
B, A = signal.butter(N, Wn, output='ba')
df_filtered["ROLL"] = signal.filtfilt(B,A, df_filtered["ROLL"])

x_axis='Date'
plt.rcParams["figure.figsize"] = (17,10)
fig1, (ax1,ax2,ax3) = plt.subplots(nrows=3, ncols=1)

xlim_start=min(df["Date"].iloc[0],df["Date"].iloc[0])
xlim_end=max(df["Date"].iloc[-1],df["Date"].iloc[-1])
# xlim_start=datetime.datetime(2018,7,3,12, 51)
# xlim_end=datetime.datetime(2018,7,3,12, 59)
xlim=[xlim_start,xlim_end]

ax1.plot(df[x_axis],df['ROLL'], 'r.')
ax1.set_ylabel("ROLL")
ax2.plot(df_with_outliers_trimmed[x_axis],df_with_outliers_trimmed['ROLL'], 'r.')
ax2.set_ylabel("ROLL with outliers trimmed")
ax3.plot(df_filtered[x_axis],df_filtered['ROLL'], 'r.')
ax3.set_ylabel("ROLL filtered")
ax1.set_xlim(xlim)
ax2.set_xlim(xlim)
ax3.set_xlim(xlim)

  csv file:
  https://file.io/yqf3ei


texts: Disparity in phase calculation between STFT and Hilbert transform

texts:
I have a signal which includes several frequencies (EEG), and I want to calculate the phase of this signal in a particular frequency in a particular moment.

For this I have used two approaches a Short Time Fourier Transform (STFT) and a Hilbert transform:

I'm not worried about a possible constant (or close to constant) phase difference between the two, but about a non constant difference, which is the case, the two phases do not correlate very much.

So my impression is that both differ for the reasons below with a longer explanation on the method:

My signal has a frequency sampling of 2ms.

Let's call the time of interest T0. And foi our frequency of interest (10 hz for this particular case, though it changes).

For a particular signal segment (800 ms) I calculate an spectrogram and find the closest peak in the spectrogram to the foi, and then I calculate the a small bandwidth around it let's call it Hilbert Bandwidth (~+- 1.25Hz). (I do this to make sure there is some signal to filter later).

For the STFT,  and I use a 200 ms window centered over T0. Therefore 100 points (given 2 ms sampling) . Giving me a frequency resolution of 5 Hz. As I want the phase in a particular moment, I tried to use a window as small as possible. Then I calculate the of the foi based ont his transform.

For the Hilbert Transform I apply a bandpass Butterworth filter (order 2) using the Hilbert Bandwidth, and use a Hilbert transformation to calculate the phase.

This two phases (Hilbert and STFT) differ, so my reasoning here is:
- That the STFT phase for the foi has contributions of closer frequencies that I cannot separate due to the length of the window. While the Hilbert does not as I have applied a bandpass filter narrower than the resolutuion of the STFT.

So I have three questions:

Is filtering in a narrow band and calculating the Hilbert transform a senseful way of calculating the phase of the signal?
Does my explanation of the difference between the two methods make sense?
For a given signal,  which of the two methods would be more reliable for estimating the phase at a particular time point in a particular frequency?


texts: What is meant by "Group delay"?in simple words?

texts:
Although likes of this question have been asked many times on DSP SE

But i am unable to understand those and i wish to have a crystal clear explanation in simple words with example

The only thing that i am able to understand from those questions is that group delay is negative derivative of phase with respect to frequency. But what does negative derivative means here and why not positive derivative here?


texts: What is meant by "information encoded in time domain"?

texts:
Author on pg 292 claims the windowed-sinc "is not a filter for signals with information encoded in time domain", because its "step response has overshoot and ringing". That phrase, along "encoded in frequency domain", are used in text, but I'm unsure they're ever explicitly defined. My guess:

Frequency-domain encoded: periodic processes in a system - something that recurs, persists the duration of the signal, and is sensibly described by Fourier coefficients
Time-domain encoded: the opposite - time-local events in a system, such as a sharp rise, temporary change in DC offset, momentary noise injection.

A single square pulse would qualify for latter, as the physical system that generated it almost certainly didn't synch hundreds of sinusoidal sources just right. A noiseless EEG would qualfiy for former, as neuronal excitations recorded as polarity alignments have periodic traits.
Is this accurate? If so, how is a filter's step response relevant in describing time-local events - and are there other "responses" characterizing time-encoded information?

Note: I'm aware of the step response's role in control systems, circuit design; I ask specifically about non-periodic "events" in a signal (if that's what time-encoded means).


texts: How do I time-stretch a signal from the frequency domain?

texts:
I have $FFT(signal)$, and want to find $operation$ such that $$IFFT(operation(FFT(signal),2)) = timestretch(signal, 2)$$
Where $timestretch$ stretches the signal in the time domain - maintaining the frequencies, unlike basic subsampling stretching which leads pitch change.
One option I know I could do is convert to a 2D spectrogram and stretch this "image" and then reverse the spectrogram, but this seems both inaccurate and expensive.


texts:
The Nyquist theorem states that the sampling rate must be twice the highest frequency to be observed. How does the length of the interval play into this relationship?
The main resource that I found here states that

For a waveform to be fully characterized, it must complete at least one full cycle within the duration of the time record.

Is this the only restriction, assuming that my sampling rate exceeds the Nyquist rate? For example, assume that the frequency band I wish to measure is from 8Hz to 13Hz. I have a sampling rate of 2000 Hz, which is way beyond the Nyquist frequency. As long as my time window is at least 1000ms / 8 Hz = 125 ms, should I be able to reconstruct the signal? For context, I am planning on running an FIR filter followed by a Hilbert transform on these windows to extract phase information, so I want the analytic signal to be mathematically sound.


texts: How can I get the continuous-time transfer function coefficients (or poles and zeros) from the corresponding discrete-time TF and vice versa?

texts:
Let's say I have a continuous time transfer function which I know its numerator coefficients $(B^c = [b^c_m, ..., b^c_1, b^c_0])$ and denominator coefficients $(A^c = [a^c_m, ..., a^c_1, a^c_0])$. Also I know its zeros $(Z^c = [z^c_m, ..., z^c_1])$, and poles $(P^c = [p^c_m, ..., p^c_1])$ and gain $K^c$.

How can I get the numerator $(B^d = [b^d_m, ..., b^d_1, b^d_0])$, denominator $(A^d = [a^d_m, ..., a^d_1, a^d_0])$, zeros $(Z^d = [z^d_m, ..., z^d_1])$, poles $(P^d = [p^d_m, ..., p^d_1])$, and gain $K^d$ of an equivalent discrete-time transfer function (considering sampling time of $T_s$)

I know about the transformation methods (like Tustin or bilinear approximation), but those transformation methods convert the transfer function from continuous to discrete and it is not obvious from that how the numerator and denominator coefficients change. In other words I need the implementation of two functions like below:

$$[A^d, B^d] = \text{function_1} \,(A^c, B^c)$$
$$[Z^d, P^d, K^d] = \text{function_2} \,(Z^c, P^c, K^c)$$

I know for function_1 there is a Matlab function called c2d() which does that, but that is not compatible with Matlab coder, and I need to write my own function. So that is why I am interested to know the explicit mathematical formulas that does the conversion between coefficients (or zeros and poles) and not the transfer function.

For the function_2 maybe I have the answer but I am not sure it can handle all cases (including single/multiple poles or zeros):

$Z^d = e^{Z^c \times T_s}$

$P^d = e^{P^c \times T_s}$

$K^d = 1-e^{-K^c \times T_s}$    ---> (is the minus sign right?)


texts: Is there a linear filter whose step response is an S-curve?

texts:
Is there a linear filter whose step-response (or ramp-response) is an S-curve?

I think the properties I'm looking for here are (1) time-reversal symmetry, and (2) at least a continuous first derivative.


texts: Second moment ergodicity of gaussian random process

texts:
How can I prove that a WSS Gaussian stochastic process with mean 0 is
mean-square ergodic in the second moment if and only if:

$$\lim_{n \to \infty} \frac{1}{n}\sum_{k=0}^n r_{xx}^2(k) = 0$$

When $r_{xx}(k)=E[x(n)x(n+k)]$ and $x(n)$ is the process.


texts: 2-d circularly symmetric low-pass filter

texts:
For a square pixel grid, the ideal 2-d low-pass filter with a horizontal and a vertical cut-off angular frequency $\omega_c$ in radians has an impulse response (kernel) $h_{\small\square}(x, y)$ that is the product of a horizontal and a vertical stretched and scaled sinc function, with $x$ and $y$ the integer horizontal and vertical pixel coordinates:

$$h_{\small\square}[x, y] = \frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c x}{\pi}\right)}{\pi}\times\frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c y}{\pi}\right)}{\pi} = \begin{cases}\frac{\omega_c^2}{\pi^2}&\text{if }x = y = 0,\\\frac{\sin(\omega_c x)\sin(\omega_c y)}{\pi^2 x y}&\text{otherwise.}\end{cases}\tag{1}$$

If $\omega_c = \pi$, the kernel of Eq. 1 is simply:

$$h_{\small\square} = [1]\quad \text{if}\quad\omega_c = \pi.\tag{2}$$

For real-valued $x$ and $y$, Eq. 1 would not describe a circularly symmetric kernel, and neither is the square-like frequency response $H_{\small\square}$ circularly symmetric. Sometimes an isotropic kernel may be desired, for example if we want to remove the effect of the choice of the direction of the image axes.

What is a kernel $h_\circ[x, y]$ such that the frequency response $H_{\circ}$ is $1$ inside a circle of radius $\omega_c$ (angular frequency in radians) and $0$ outside it?


texts: Simple way to replicate bloom (glow) effect in image using Matlab?

texts:
I have some simulated images of light bulbs in Matlab. What I want to do is replicate the bloom (glow) effect on these images. The degree of the effect (i.e. how much the light bulbs bloom) should be controllable.

Is there any implemented function in Matlab Image Processing Toolbox help me to achieve this. Or is there any simple algorithm for this. I don't want a sophisticated algorithm since I don't require the bloom effect to be precisely replicated.

I would appreciate any suggestion.

Thank you.


texts: Which image transformations are considered continuous transformations?

texts:
In the article J. Matas, O. Chum, M. Urban, T. Pajdla: Robust Wide Baseline Stereo from Maximally Stable Extremal Regions it is stated that:

  The set (of extremal regions) is closed under continuous (and thus perspective) transformation of image coordinates, ...

I know what perspective transformations (and consequently, affine transformations) are, but could anyone give examples of other transformations that fall under the category of continuous transformations of image coordinates? How would such transformations be characterized?


texts: What is the purpose of wrapping the negative times of a response function in discrete convolution?

texts:
I am trying to rationalize a figure given in the Numerical Recipes in C in the section of Fourier based convolution and deconvolution.
The authors show the example of a continuous convolution with a zero centered response function r(t). Convoluting the continuous signal blurs out the edges and there is no shift.
What is somewhat confusing is that they suggestion for the discrete convolution analogue of the continuous convolution, we should wrap the negative time section of the response function to the extreme right end of the array. What is the mathematical justification of doing N/2 circular shift, where N is even number of samples. The authors do not mention as N/2 circular shift, but it seems this is what they are doing.
These are their actual figures.

It seems the authors are suggesting circular shifting the zero centered response by N/2 (N is even) to get an equivalent convolution to continuous convolution. Why can't we have discrete version of the response function as shown by the squared markers which I drew in MATLAB? I do not get equivalent convolution results if we convolute a square wave with the sampled response function on the right and the left as suggested in Numerical Recipes.


texts:
I'm not quite sure if I can ask about a vendor-specific hardware here, please teach me if not.

We have a TI's DSP EVM, one of EVMK2Ls (XEVMK2LX/XTCIEVMK2LX) and wanna do some digital signal processing using it.

To get analog signal data, we purchased a 3rd-party ADC EVM, DEV-ADC34J22, considering necessary specs such as sampling rate and support for the high speed transmission protocol, JESD204B.

The ADC EVM, however, doesn't have a connector for direct connection to the DSP EVM (which has FMC connectors), instead has the HSMC connector only, so we had to buy an adapter card, KY-FMC2HSMC.

Now, we can connect the two EVMs, DSP-EVM and ADC EVM, to each other, but don't know how to use the combination for our purpose. As mentioned above, we want to get data from the ADC and do some processing on them in a programmatic way (in C, more specifically). We posted many questions on the TI's E2E community forum and got some answers, and finally concluded that we can't achieve the goal with the hardware configuration (even if it's theoretically possible, practically not possible for us).

So, now we are considering to purchase another ADC EVM, which seems to be TI's 1st party EVM and to connects to the DSP EVM directly with the FMC connector. It has the same chips (ADC/LMK) made by TI with the ADC EVM we bought, and therefore meets our requirements. We hope we can just use it out of the box.

We got some advices from TI's engineers at the community forum actually, but the options suggested are not affordable for us due to the limited budget. We asked about the TI's ADC EVM we wanna try, but didn't get the answer yet. We came to know this DSP stack exchange site today and would like to get some advice also from experts who is not a TI employee, even though the answers from the TI's engineers were also helpful.

Any other recommendation on the hardware configuration would also be appreciated. We just want to be able to acquire and process a large number of data programmatically as quickly as possible. The number is upto 20M, and we need to do FFT and other transforms on the data using a DSP or anything having enough processing power for that.

Thank you very much!


texts: What does the sample entropy of a digital signal mean?

texts: What sample entropy signifies about a signal and how it used to analyse a signal?


texts:
I am using matched filters using GNU Radio to extract from noise a CW Morse signal. The filter actually performs pretty well in various SNR environments. My problem is that the threshold upon I decide for a presence or absence of a CW transmission vary with different SNR/CW frequency/WPM/incoming signal amplitude setups.
Conventional normalization methods, (eg dividing with the max value) are not applicable since the program runs continuously and the max value is not known a priori.
I believe that a moving normalization scheme should deal with this problem. In the past, for preamble detection I have used the moving normalized autocorrelation method that is depicted below (C and P are moving sums). I tried to adapt this approach in my case but with no luck.
My question is, if there is a method to normalize either the output signal power of the matched filter or the input signal itself, so I can rely on a one predefined decision threshold for CW signal presence/absence.


texts:
I saw a couple of definitions for minimum-phase in different textbooks and I'm trying to understand what the implication of each of them.
The first definition I saw was: An invertible system which both it and it's inverse are causal and (BIBO) stable.
The second definition is: An invertible system which both it and it's inverse are causal and have finite energy.
It was written in one textbook that had the second definition as an alternative to the first one that these constraints allow the system to have zeros and poles on the unit circle and I don't understand why (as opposed to the first def).
just to make sure, I'm not restricting myself to rational systems.
Another thing I saw in some textbook and didn't understand is that the second definition implied that the transfer function of the system and it's inverse are analytic in the exterior of the unit circle.
So, to conclude, my questions are:

why we have both of the definitions of minimum-phase?
why the 2nd def implies that the system can have zeros and poles on the unit circle as opposed to the first def?
why the 2nd def implies that the system and it's inverse have analytic transfer functions in the exterior of the unit circle?

Sorry for any grammar mistakes and thanks for any clarifications.


texts: Digital Lowpass butterworth filter with cut off 500Hz and sampling rate 1.25MSPS

texts:
I am trying to simulate a distributed sensing system and I need to filter only frequencies lesser than 500 Hz(low pass filter) from acquired signal with sample rate 1250000 Samples/sec using below mentioned program:

 %Time specifications:
   Fs = 1250000;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 0.25;             % seconds
   t = (dt:dt:StopTime-dt);     % seconds
   %Input Sine wave:
   Fc = 300;                     % hertz
   dataIn = cos(2*pi*Fc*t);

   % Plot the signal versus time:
   figure;
   plot(t,dataIn);
   xlabel('time (in seconds)');
   title('Input signal');
   zoom xon;

   wc=500; % cut off frequency
   Wn=wc/(Fs/2); % normalized frequency
   [b,a] = butter(6,Wn,'low');
   figure;
   freqz(b,a)

%Plot filtered output
    dataOut = filter(b,a,dataIn);
    figure;
    plot(t,dataOut)
    title('dataout')
    zoom xon;

However I am unable to filter the desired frequencies below 500 Hz, as the cut off frequency is dependant on the sampling rate. Kindly suggest what to be done as I am new in MATLAB.

Thanks in advance.


texts: MATLAB: Normalize Image in Frequency Domain

texts:
I have 3 different images and I have the data which is in frequency domain. The data size is : 192x256x88x32x3.

First 3 sizes are frequency domain dimensions
    4.dimension: channels
    5.dimension: different images' data in frequency domain

When I take FFT3D for this data and sum of squares in the channel direction I obtain the image. I want to normalize all this images such that their maximum image domain intensity value should be 1. The process I'm doing is taking 3D FFT for 32 channels. Then I'm finding maximum of the image intensity and multiplying the maximum value with frequency domain.

Is my approach true?
After taking absolute value, there are a
    lot of outliers. Will it be problem?


texts: Impulse response of system with $\max\left\{\lvert x[n] \rvert\right\}\ge\max\left\{\lvert y[n]\rvert \right\}$

texts:
Find a relation/condition for the impulse response, so for any $x[n]$ (input) with $y[n]$ (output) this relation is true:
$$ \max\left\{\lvert x[n] \rvert\right\} \ge \max\left\{\lvert y[n]\rvert \right\}$$
The system is discrete, linear and time invariant.

I tried to replace the $y[n]$ with the convolution formula,but it doesn't work. I tried to use Fourier transformation, but it doesn't work either. Any suggestions?


texts: Using the linear operator to check for time invariance of a differential equation?

texts:
I have a differential equation $$\frac{d^2y(t)}{dt^2}+y(t) = \frac{dx(t)}{dt} + x^2(t)$$ and I need to see if this system with input $x(t)$ and output $y(t)$ is time invariant and linear. I tried to break up the equation by substituting a linear operator D as the derivative.

So I get: $$y(t) = \frac{Dx(t)+x^2(t)}{D^2 + 1}$$

But the problem is, this equation is meaningless since you are dividing by the derivative operator. However, mathematically, it works out, so can I get $y(t)$ in this form and test for linearity or time invariance using this method?


texts: Calculate the time delay introduced by group delay for IIR-Filters

texts:
First of all, let me show what I mean by a simple example of an FIR-Filter:

I designed an FIR-Filter with a constant group delay of 323 samples and my sampling rate is 500 Hz. Now I can calculate the time delay by this filter with 323/(500 Hz) = 0.646s.

An IIR notch filter has a frequency dependent group delay and I am interested in the time delay around 50 Hz.

We know the group delay is around 2.5 samples at 50 Hz and we sample with 500 Hz. Can I do the same calculation as before? This means the time delay for 50 Hz is 2.5 samples/500 Hz = 5ms.

I am not sure but I think the correct calculation is with fs/2, that means it's 2.5/250 = 10ms.


texts:
A signal $x[n]$ is to be transmitted over a communication channel. The communication channel is described by an FIR filter of length $M$ such that the received signal
  is given by
  $$r[n]=\sum_{k=0}^{M-1} x[k]h[n-k]\text.$$

  The data are transmitted in blocks where each block contains N samples of a data
  stream $d[n]$, and in addition each transmitted block also contains a so-called cyclic
  prefix where the last values of the block are duplicated at the beginning of the
  block. Consider for example the first block of data ${d[0], d[1], . . . , d[N − 1]}$.

  The
  corresponding block to be transmitted with the cyclic prefix is
  ${x[n]}_{k=0}^{N+K-1}= {d[N − K], d[N − K + 1], . . . , d[N − 1], d[0], d[1], . . . , d[N − 1]}$
  $i.e. x[0] = d[N − K], x[1] = d[N − K + 1], . . ..$
  At the receiver end, the first $K$ values are discarded and the next $N$ values are saved in a buffer  ${y[n]}_{k=0}^{N+K-1}={r[K], r[K + 1], . . . , r[K + N − 1]}, i.e. y[0] = r[K], y[1] =r[K + 1], . . ..$
  a) Show that at the receiver end the original data stream d[n] can be recovered
  using the following algorithm if the cyclic prefix is sufficiently long.

  $Y [k] = FFT({y[n]})$
  $D'[k] = Y [k]/H[k]$
  $d'[n] = IFFT (D'[k])$

  Show all workings, and in particular determine the length $K$ of the cyclic prefix such
  that $d'[n] = d[n]$,$ n = 0, . . . , N − 1,$ and show exactly how $H[k]$ should be computed
  from $h[n]$.

To solve this question for the given info do we have to convert the sequence into odd and even sequence respectively and then arrive at the FFT or is there any alternative? I don't need the solution. A hint would be good to start.


texts: Beginner level: deconvolving Gaussian noise to find initial PDF

texts:
noob question.

An initial process has an unknown pdf, P(X), which is then subjected to additive Gaussian noise.  Repeated sampling of the sum is performed so that

$y=x+noise$

is repeatedly measured and the pdf for y is determined, P(Y).  Is there anything that can be said analytically about P(X)?

I'd appreciate any simple references.

Thank you


texts: Stabilizing the inverse transform of a system

texts:
Currently I am trying to understand an unknown mechatronic system. I used white noise as an input signal and measured the output of the system. Then I estimated a second order SISO transfer function (Output/Input) in the frequency domain using Matlab.
sys=Out/In
h = idfrd(sys, freq*2*pi, Ts);
np = 2;
nz = 2;
Gs = tfest(h, np, nz)
opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', 2*pi*120);
Gz = c2d(Gs, 1/4000, opts);

$$
Gs=\frac{3.765\cdot s^2 - 528.7 s - 5.113 \cdot 10^4}{s^2+183 \cdot s + 4.117 \cdot 10^5}
$$
A bode Plot of the derived second order system $Gs$ (blue) and the measured data (orange) shows good agreement in amplitude and phase above $60\;Hz$.

Further, a pole-zero plot of $Gz$ shows that all poles are stable (inside unit circle of z-plane), but that there is one zero outside of the unit circle of the z-plane.

Now, I want to counteract the behaviour of the mechatronic system in a way that I achieve a unity output. So by modifying the input to the system I want to basically get rid of the system characteristics and achieve a pure unit transfer behaviour.
IGz=1/Gz; % so that IGz * Gz = 1

Mathematically, I need to multiply the input signal with the inverse transform of the system. Here a bode plot of the system $Gs$ (blue) and its inverse transform $IGz$ (orange) basically shows what I want to do.

But because of the zero outside the unit circle, the inverse transfer function now has a pole outside the unit circle and thus is unstable. Therefore I wonder how I can counteract the behaviour of the system while not getting an unstable system?

Attempt
So far, I have tried to replace the unstable pole by reflecting the pole inside the unit circle as advised here: $(z-p)$ by $(z-1/p)$ (as the unstable pole is purely real I don't need to do the complex conjugate). This keeps the amplitude response, but changes the phase response (blue: unstable $IGz$, orange: stabilized $IGz$):

However, I would like to keep the phase response as close as possible to my inverse transfer function as well. But, as I don't care so much about what happens below $60\; Hz$ I wonder if there is an approximation of my inverse system which is good above $~60\;Hz$ in amplitude and in phase, but can deviate below $~60\;Hz$ and is stable. I wonder if there is an established way how one would derive such a system?
Here is what I have tried so far.
pIG=pole(IGz);
zIG=zero(IGz);

z=tf('z',1/4000);
Test=((z-zIG(1))*(z-zIG(2)))/((z-1/pIG(1))*(z-pIG(2)));
Fac=evalfr(IGz,100)/evalfr(Test,100);
Test1=Fac*Test;
bode(IGz,Test1)

The result looks good at higher frequencies, however the phase could be better below $100\;Hz$ or so (blue: unstable $IGz$, orange: stabilized $IGz$ with phase kept roughly the same above $100\;Hz$):

Edit
I just looked at the impulse responses. The system $Gs$ of course is oscillating, but decaying.

However, the impulse response of the inverse system is, of course, unstable.

The impulse Response of the second stabilized system is still overshooting with high amplitude. I wonder if this works, since I, of course, am limited in output w.r.t. the controller.

With the input $u$ to the controller and the output $y$, I want to do the following $u \rightarrow IGz (controller) \rightarrow y \rightarrow Gz (mechatronic\; system) \rightarrow y_2$. During all this, the output $y$ of the controller can't get too big and it should obey $\frac{y_2}{u}=1$.
How can I best come up with an $IGz$ that achieves this?


texts: What should be the singal duration in time doamin to correctly measure its frequency using FFT

texts:
I am generating signal consisting of 4 frequencies 380,460,500,580 Hz using 20500 as Fs. I generate this signal for different durations of 0.004s,0.008s,0.02s,0.04s.

Fs=20500;
values=0:1/Fs:duration;
sig1=sin(2*pi*FREQ*values); % FREQ can be 380,460,500,580
sig2 , sig3, sig4
sigt=normr(sig1+sig2+sig3+sig4);

When I take FFT I get different results for different time durations.

   NFFT=2^nextpow2(length(sigt));
   Y=fft(sigt,NFFT)/length(sigt);
   f=Fs/2*linspace(0,1,NFFT/2+1);
   figure, plot(f,2*abs(Y(1:NFFT/2+1)),'r')

for 0.004s, I get one peak at 480.5 Hz
for 0.008s , I get one peak    at 480.5 and two lower side peak at
320.3 and 640.6.
for 0.02s ,I get    four peaks at 360.4,440.4,520.5,600.6
for 0.04s , I get four Peaks at 380.4, 460.4,500.5,580.6

All frequencies have time period less than any of above durations. Is there any specific criteria that signal must be in some length so that FFT can detect its frequency.
Is there any other method exist which can measure frequency for such short durations?
Please share any knowledge.


texts: How to filter uniform noise in a microcontroller?

texts:
By uniform noise I mean a band of a fixed width, no spikes, etc.  This has to be one of the simplest examples. I'm looking for the appropriate theoretical approach as much as something relatively computationally efficient.

The true signal changes slowly (slower than 1Hz, a battery voltage).  The noise band remains fairly constant (~50mV noise on a ~4V signal).

Initial ideas:

Long moving average (e.g. a 32-point moving average).
Leaky-integrator/low-pass filter/FIR filter:
"Low-pass filter" in non-EE, software API contexts
Recursive arithmetic average:

y_new = (new_sample + y_last) / 2

I'd love to learn more about what I think is a geometric average.  I'm not sure I can quantify the time-to-settle/time-to-seed for these options either.

If there are more standard approaches, I'd love to learn about those as well.


texts: So if I were to have an $n$ length sample of a given sample rate and were to run it through an FFT, resulting in $n$ complex elements. How would I obtain the frequency bins, such that if I know the frequency of some noise, I can remove the noise?


texts: SNR in case of complex channel and fading

texts:
I am trying to estimate SNR mathematically given in a paper.

$$ y = hx + n$$
$h$ is channel fading, that is Rayleigh in my case, while $n$. Both $h$ and $n$ are random variables.

 n = a * (randn(1,1) + 1j*randn(1,1));        %complex gaussian
 h = abs(b * (randn(1, 1) + 1j * randn(1,1)); %complex rayleigh

Here in this case when both $h$ and $n$ are complex, is $a/b$ is my SNR?

is $a$ the variance of noise?
is $b$ the channel gain?


texts: DFT of complex signal using Goertzel algorithm in C

texts:
I am implementing a BFSK frequency hopping system with TX and Rx modules. I am using Goertzel Algorithm at the receiver end to demodulate the data i.e. to determine the carrier frequency of the recevied signal data. Following is the implementation of the goertzel implementation:

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* modData)
{
int k,i;
float   floatnumSamples;
float   omega,sine,cosine,coeff,q0,q1,q2,result;

floatnumSamples = (float) numSamples;
k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * 3.1416 * k) / floatnumSamples;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

q0=0;
q1=0;
q2=0;
result=0;

for(i=0; i<numSamples; i++)
{
    q0 = coeff * q1 - q2 + modData[i];
    q2 = q1;
    q1 = q0;
}

result = sqrtf(q1 * q1 + q2 * q2 - q1 * q2 * coeff);
return result;
}

The implementation works well for a signal where the received signal is a real signal ( a set of float values). Would this algorithm still work for complex signals. The complex signal would be typically be represented by a set of float values, with each alternate element representing a real and imaginary value. In this case, how should the above goertzel function be used, should it be used to calculate the output for the real and imaginary components separately, or should it read the entire complex signal array elements in a sequence (both real and imaginary array elements) and then  calculate the magnitude at the end of it. If not, can some please point me to how to get the goertzel output for a complex signal samples. I tried googling it but could not find anything concrete about how to use goertzel algorithm on complex signals.

Also, I wanted to ask one more thing. If my signal has carrier frequencies at both 4000 and -4000 Hz and a sampling frequency of 18000 Hz, then for a real signal as input, the goertzel output is the same for input frequency of 4000 Hz and 14000Hz (this is what should be). But since 4000 Hz and -4000Hz are one of my carrier frequencies, I am not able to determine whether the received signal is at 4000 Hz or -4000 Hz. Is there any method to do this without doing a frequency shift. I am able to differentiate using frequency shifting (multiplying the received signal with a sine wave before demodulating). But I am looking for a method to do this without using the frequency shift technique.

--Anshu


texts:
One of the approaches I have been searching so far is the application of utilizing Neural Network to classify Spectrogram images, which is heavily used in speech recognition.

I am conducting an experiment to identify the similarities of binary executable (raw 0s and 1s), and planning to deploy the same concept of speech recognition on binary executable. I have been preparing the executable with the following criteria:

1- The sampling rate is 8 KHz.

2- The bit Depth is 8-Bits. This means that each sample's amplitude will be represented using 1 Byte.

I am being stuck in generating the Spectrogram. Now my background is networks and programming and I have been doing a lot of research to generate it. The challenge is the identification of the frequencies of the samples I have in order to have the 3-dimensions of spectrogram (frequency, time and amplitude). I did a research, where FFT can be used to generate the frequency domain from the time domain, but not sure if it is feasible for binary executable.

Python supports .wave files spectrogram plotting, but the files I have are executable and not music or sound.

1- Can I identify the used frequencies in my binary by having the sampling rate and amplitude?

2- Can I generate a spectrogram of an executable binary, even if it is NOT originally speech or music?

3- Is there any tool/code that I can utilize for the spectrogram generation?


texts: What assumptions should be used to invert spectrum into time domain data?

texts:
This question comes from a realistic problem that how to obtain a time domain wind pressure from known wind spectrum. Or similarly, gets rail surface curve from spatial rail spectrum.

Obviously, this question does not have unique solution, since phase information is not given. Usually we assume that phase of the time/spatial domain data in uniformly distributed in ($-\pi$, $\pi$], and using a random distributed phase to do inverse Fourier transform to obtain the time/spatial domain data. However, I feel not so comfortable to assign such rapid changing phase to the spectrum. I wonder whether there are some more reasonable assumptions to get these jobs done?

BTW, question ifft of signal without phase information discussed some strategies on how to do this job. @LutzL suggested Human perception can be used to help determine the phase information. But I cannot understand how Human perception is working. @Hilmar said minimum phase may be a choice, but I do not know whether there is a physical meaning of minimum phase. Is there any discussion to show minimum phase is a good choice for physical phenomena, like wind.


texts:
I am applying a Savitzky-Golay filter to a signal, using the scipy function.

I need to calculate the lag of the filtered signal, and how much is it behind the original signal.

Could someone shed some light on this matter? How could I calculate it with scipy? How should I interpret the result correctly?

I would be very grateful!


texts: I'm new to OFDM modulation. During OFDM modulation in the frequency domain, are the negative frequencies modulated with a conjugate version of the data of the positive frequencies, or do the negative and positive frequencies carry completely different data?


texts: Detecting frequency when frequency is not sufficiently greater than baud rate

texts:
I'm writing software this software for decoding audio samples of FSK data back to bits (ASCII, etc). The software is working really well at this point for decoding FSK when the mark and space frequencies are sufficiently greater than the baud rate.
However, I detect a frequency of $0\textrm{ Hz}$ when the mark and space frequencies are decreased to values which are not sufficiently greater than the baud rate. I've been using a pretty basic zero crossing algorithm for detecting the frequency present in a given interval. I gather that this is because the frequency is not great enough to have any zero crossings during the interval that I sample.

What are some other methods for detecting whether the mark or space frequency are present in a given interval when the baud rate is not sufficiently greater than said frequencies?

And as a follow-on question, is it common for the mark and space frequencies to not be sufficiently greater than the baud rate in real-life FSK implementations?
For example, I read (here and in other places) that the mark and space frequencies for Bell 202 (1,200 baud) are $1200\textrm{ Hz}$ and and $2200\textrm{ Hz}$, respectively -- the algorithm I'm currently using can detect the $2200\textrm{ Hz}$ frequency but not the $1200\textrm{ Hz}$ frequency.


texts: Difference between "scaled white noise " vs. "unscaled white noise"

texts: I have been searching Matlab user-guide for Wavelets. In their graphical user interface for denoising, there is an option for "unscaled white noise" and "scaled white noise" for the noise structure. What is the main difference between the two? Does the scaled white noise imply that the noise is amplitude changing with time? I have searched Google books and Google Scholar, nobody defines it. Thanks.


texts: Time-domain LPF not showing expected behavior

texts:
I am trying to implement a simple first-order Butterworth Low-Pass filter in Python. I have some code that makes use of scipy.signal.butter and scipy.signal.filtfilt. It works fine, but I wanted to make sure I could build the filter myself as well.

I start with a bilinear transform $s = \frac{2}{T_{s}}\frac{z-1}{z+1}$ on my first order LPF $\frac{\omega_{c}}{\omega_{c}+s}$. Making that substitution and converting from the $z$ domain to the time domain gives me the equation that I'm trying to reproduce:

$y_{n} = \frac{1}{\omega_{c}T_{s}+2}[\omega_{c}T_{s}(x_{n}+x_{n-1})+(2-\omega_{c}T_{s})y_{n-1}]$

I'm confident that this expression is correct, but my attempt at coding it doesn't seem to be working.

EDIT: I changed my implementation so what's actually happening is more clear. As a side effect, I'm now seeing a different type of strange behavior than I was previously. The new code is below:

T = 2                                                  # time length over which samples are collected (s)
f0 = 5                                                 # fundamental frequency (Hz)
fs = 1000                                              # sampling frequency
Ts = 1/fs                                              # sampling period
fc = 10                                                # filter cutoff frequency
wc = 2*pi*fc                                           # cutoff frequency in radians

# Create wave input to be sent to filter

t = np.linspace(0,T,T*fs+1)
x_storage = [10 + cos(2*pi*f0*x) for x in t]
plt.plot(t,x_storage)

# Create a place to store output data

y_storage = np.zeros(len(t))

# Initialize input and output data

y = 0                                                  # y_(n)
yprev = 0                                              # y_(n-1)
x = 0                                                  # x_(n)
xprev = 0                                              # x_(n-1)

# Begin calculation

for idx in range(len(t)):
    # One timestep of the calculation
    xprev = x                                          # move the old x_(n) value to x_(n-1)
    x = x_storage[idx]                                 # bring in the current datapoint to be the new x_(n)
    yprev = y                                          # move the old y_(n) value to y_(n-1)
    y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts)    # calculate the new y_(n) according to the filter equation
    y_storage[idx] = y                                 # store the y_(n) value for later plotting

 plt.figure()
 plt.plot(t,y_storage)

My input is a wave with a DC offset and a frequency component below the cutoff threshold, so I would expect the output wave to look more-or-less the same. Instead, I see the DC bias vanish and the signal itself get hugely attenuated (see below). Any thoughts would be much appreciated!


texts: When simulating a voltage divider using wave digital filters I get inversed polarity outputs

texts:
I'm studying wave digital filters for my bachelor's thesis and I'm finding problems in implementing a voltage divider to test my knowledge.

Here's the circuit that I'm trying to simulate and the WDF tree that I think it's valid

To simulate the circuit above I created the following matlab code

clear all

t = [0:0.1:10];
f = 10000;
fs = 100000;
n = length(t);

% Resistance of resistors (duh)
Rv = 10e3;
R1 = 10e3;
R2 = 10e3;

% Voltage input
V = sin(2*pi*f/fs*t);

% Initialize arrays
Ac = zeros(1,n);
As2 = zeros(1,n);
Br2 = zeros(1,n);
Ar2 = zeros(1,n);
Bs2 = zeros(1,n);

% Series connector port resistances
Rps1 = Rv;
Rps2 = R1;
Rps3 = R2;

% Series connector port scattering parameters
lps11 = (2*Rps1) / (Rps1 + Rps2 + Rps3);
lps12 = (2*Rps2) / (Rps1 + Rps2 + Rps3);
lps13 = (2*Rps3) / (Rps1 + Rps2 + Rps3);

% Iterates over the whole input wave
for i=2:n
%     Inputs waves to the series connector
    As1 = V(i);
    As2(i) = 0;
    As3(i) = 0;

%     Reflectes waves to the series connector
    Bs1 = As1 - lps11*(As1 + As2(i) + As3(i));
    Bs2(i) = As2(i) - lps11*(As1 + As2(i) + As3(i));
    Bs3(i) = As3(i) - lps11*(As1 + As2(i) + As3(i));
end

% Calculates output voltages over resistors
Vr1 = (Bs2+As2)/2;
Vr2 = (Bs3+As3)/2;

% Plot outputs
subplot(2,1,1);
plot(t,V);
subplot(2,1,2);
plot(t,Vr1);

Which gave me the following output

The first plot is the input voltage and the second is the output voltage over one of the resistors

As can be seen the magnitude of the output is the same as expected (1/3 of the input), but the signal is multiplied by -1.

Is there something fundamentally wrong that I'm doing so that even the simplest circuit can't be correctly simulated?


texts: FFT response for noisy signal, averaging the response and

texts:
I try to analyze a system's state with acceleration sensors by measuring the frequencies and amplitudes.

Since a cheap sensor is used, I have a noisy signal in the beginning (up to 40Hz - see below). That is the first thing that I don't get. With another sensor (much higher resolution but also different position) this noise couldn't be observed.

Do you have any idea where this might come from?

Due to the noise and to reduce uncertainty, I perform multiple measurements which I intend to average afterwards. The second question would, therefore, be:

Where do I best average the measurements in order to obtain the "correct" amplitude/frequency? I have the following ideas:

average the signals and perform a single fft. The maxima result from the            maximum of the transformed signal
average the signals after the transform and take the maxima of the resulting signal
average only the maximum/frequency of the single transforms of each measurement

Furthermore, I observe a sort of frequency shift within the measurements (the left vertical line is the excitation frequency, while the right one is the response of the averaged fourier transoforms). Here, I excite the system with 183 Hz and expect the maximum at this very frequency (or in a really small bound around it). Instead, the responses occur to be around 189 Hz and sometimes even more. Due to my lack of experience I tried to conduct some tests to "check" if my algorithm works correctly:

I generated a sine wave and added different offsets and noises. The frequencies would be met with an accuracy of 0.15 Hz (see below). Therefore, I guess it's not the algorithm.

Do you have any experience with this procedure and how much shift in frequency is normal?

Do you have any suggestions for further (pre-) processing of the signal? I was thinking about a highpass-filter to get rid of the noise in the beginning. Also, I think it might be interesting to apply a bandpassfilter only around the "bound" and analyze the signal in the time-domain - what do you think about that?
I am curious to find out if it's possible to observe the switch of a valve position. Since this is a transient behaviour I am wondering which method to apply. I'd try a PSD and to cut the evaluation of the FFT into pieces and perform it for smaller periods. Do you have any ideas?


texts: Deriving statistics of band limited Random Noise

texts:
The question is:

Consider a continuous random number with a Gaussian distribution of mean $\mu$ and variance $\sigma$ . The RV is measured from time $t=-\infty$ to $t=\infty$. This time domain signal $x(t)$ is passed through a first order low pass filter with cutoff at $f_o$. What would the statistics of the output be? Would it still be Gaussian? If yes, then what is the new mean and variance?

I am unable to understand where to start with this. Any help is greatly appreciated.


texts:
I want to test some code generated by this site. I selected Bessel LP 1st sample rate $600\textrm{ Hz}$ corner $8\textrm{ Hz}$ long $10\textrm{ bit}$. If I adjust the code for octave to be:

clear all;
close all;

N = 1024

v0=256;
v1=256;

m = [];

for mi = 1:128
    ip=round(rand(N,1)*512);
    for i = 1:N
        v0 = v1;
        tmp = ((((ip(i) * 2699550)/32)
              + ((v0 * 3856860)/2))
              + 1048576) / 2097152;
        v1 = tmp;
        op(i) = v0 + v1;
    endfor
    op = fft(op,N);
    m = [m; abs(op(2:N/2))];
endfor

x=[2:N/2];
y = mean(m);

semilogx(x,20*log(y)-170);

I get this plot:

which looks vaguely correct but the slope is $\approx 13\textrm{ dB/octave}$ and not the expected $6\textrm{ dB/octave}$.

I have limited experience with DSP so I need a method to verify the actual code. Can someone recommend a procedure for verifying code like this (that works)?


texts: If I have an image whose range is from $-18.3667\cdot10^{-5}$ to $9.3127$, how could I normalize it on the same gradient so that they fall between a specific range, say fro $0$ to $255$?


texts: Mathematical reason for MCPP in audio recognition?

texts:
MCPP uses a weird graph of human phychoacoustics measurements called MEL scale to process digital audio.

Why is MEL used in DSP voice recognition if it's based on logarithmic human perception, and why is it effective?

MEL theory gives the ratio that average humans logarithmically hear equal distance pitches, and the cause for that perception ratio is unknown.

Is it because humans speak frequencies using the same frequency ratios by which they hear them, so that PC's need to know those ratios?

Is there a physical reason for the MEL to work well in DSP? i.e. is it an unwitting compromise between equal loudness contour relative to physical sound energy and the loudest most useful frequencies in voice recognition?


texts:
If you put a wave packet through the passband of a 1st-order low-pass filter, it will be delayed by the group delay of the filter, and remain the same amplitude, right?

If you put the same wave packet through a complementary 1st-order highpass filter with the same cutoff frequency, the group delay curve is the same, so the delay of the packet will be the same, but the gain is much lower, so it will be both delayed and attenuated to negligibility.

Since the output of the highpass filter is very small, if you sum the outputs of these two filters (as in an audio crossover), I would expect it to be negligibly different from the output of the lowpass filter:  Large delayed signal + very small delayed signal = large delayed signal.

Yet if you sum the filter responses, the amplitude is 0 dB everywhere, and the phase is 0 everywhere, and therefore the group delay becomes 0, which would mean that the wave packet comes out with no delay and no changes.  I don't understand how this can be possible.  Don't filters always incur delay?  How can a filter (which also has positive group delay) undo the delay caused by the other channel, especially when this is happening in the stopband?

Which part am I misunderstanding here?

  The best-known crossover types with linear phase are first-order non-inverted crossovers, ... The first-order crossover is minimum phase when its outputs are summed
  normally; it has a flat phase plot at 0°. - The Design of
  Active Crossovers

and

  Here the result of summing the outputs together produces 0° phase
  shift, which is to say that the summed amplitude and phase shift of a
  1st-order crossover is equivalent to a piece of wire. - Linkwitz-Riley Crossovers: A Primer: 1st-Order Crossover Networks

Testing on actual pulses shows how the lowpass (blue) delays the pulse, as expected, and how the highpass (green) can combine with it to produce the original (red) pulse, but how is the highpass pulse occurring before the original if the highpass filter is causal and has positive group delay?  Intuition is failing me.

It does show that the highpass output is not as negligible as I imagined, and the delay is more negligible than I imagined, and as you move the carrier frequency around, these two properties change in a proportional way (smaller delay requires lower amplitude highpass output to correct it).  But I still don't really understand it.


texts: Distorted output of delay VST effect in C

texts:
I'm getting some odd results with some code I'm working on. Help me with the implementation of this delay effect?

The code takes an input, adds delayed input to it and then passes it to the output whilst also copying it back to the delay line. Control of the delay level and feedback are provided to control respectively the amount of delayed signal added to the output and returning to the delay line.

Here's the process of my delay effect.

void delay::processReplacing(float** inputs, float** outputs, VstInt32 sampleFrames)
{
        float *in = inputs[0];
        float *out = outputs[0];
        float delayed_sample;
        int delayed_sample_index;
        for (int frame_idx=0;frame_idx< sampleFrames;frame_idx++)
        {
                delayed_sample_index = (((int)(cBuf.index-delay_samples) + CIRCULAR_BUFFER_LENGTH) % (int)CIRCULAR_BUFFER_LENGTH);
                delayed_sample = (cBuf.buffer)[delayed_sample_index];
                //write current sample to buffer
                (cBuf.buffer)[cBuf.index] = in[frame_idx] + (delayed_sample*feedback);
                //Calculate + write output sample
                out[frame_idx] = in[frame_idx] + (delay_level*delayed_sample);
                //Increment buffer index
                cBuf.index = ((cBuf.index+1)) % CIRCULAR_BUFFER_LENGTH;
        }
}

On the whole it works pretty well and functions as a delay effect should. However, the output I'm getting is slightly distorted. I'm not sure if this is an error in the arithmetic of sample times and values or something else. Any ideas?


texts: Sine of frequency 0 contains sines of all frequencies at once in it

texts:
You know that a sine corresponds to a pulse by J.Fourier transform. The lower is the frequency, the closer is the pulse to the origin. A constant signal is a sine (or cosine, that may be important) of frequency 0. It is a pulse in the origin. This is ok, since we know very well that a pulse transforms to the white spectrum: it is a combination of all sine waves at once.

But, I see a contradiction here. I have just shown that a constant is also a sine of frequency 0 and, also, a combination of all frequencies. A sine of frequency 0 is a sum of all frequencies. How is this possible?


texts: Recursive form representation of FIR filter

texts:
Do we have methods to represent FIR filter in recursive form?

One of the methods that I know is frequency sampling method.

Can someone provide me any method?


texts: Should I sample at twice the bandwidth or twice the highest frequency?

texts:
I am confused about what rate I should sample at.  I've heard 2 different ways:
1) Sample at twice the highest frequency
2) Sample at twice the bandwidth

If I have a signal composed of just cosine, such as:
$$x(t) = \cos(4 \pi t) + \cos(2 \pi t) + \cos(1 \pi t)$$
What is the correct minimum sampling rate?


texts: Using MATLAB's Spectrogram Function for Analysis

texts:
In the documentation - sepctrogram().
Question 1)s = spectrogram(signal) and spectrogram(signal) are two commands to plot the spectrogram. However, the variable s is complex valued. I am unable to understand which output of the spectrogram is used to generate the image plot?
Question 2) How to determine the best values of the parameters window and noverlap ? Should noverlap be 50% of the signal length (number of elements in the time series) or 90% etc? What if it is zero then what does it mean? My dataset has sampling time = 1sec. I remember reading somewhere that the window should be at least roughly twice as long as the period of the lowest frequency. So, for my case is w=2 since frequency = 1?
I was thinking of using pspectrum(signal,'spectrogram') which outputs the spectrogram and use the output values as inputs to the spectrogram() function. But again, I don't know which output values from pspectrum can be used, if at all that is possible.


texts: Where can I find more information about this formula relating a window's time and frequency size?

texts:
I remember having seen a formula of the form

  $C=\frac{B}{L}$

that relates a duration $L$(s) of a window to its width $C$(Hz) in the spectrum (defined at -3dB or so). The constant $B$ would be tabulated in function of different windows.

Does anyone know how this formula is called or where I can find more information about it? I searched on Google (Scholar) but dit not find it.


texts: Jamming resistance in spread spectrum systems

texts:
In direct sequence spread spectrum, which is the basis of CDMA systems, the baseband signal is multiplied with a pseudo noise (PN) sequence which has a high chip rate, thus the resulting modulated signal will have high rate and thus larger bandwidth. But the power of the original baseband signal and the power of spreaded signal resulting from the multiplication process is the same.
Then, why this spreading process provides jamming resistance property?

Viewed in time domain, they have the same power. Does the jammer depend on the power value to detect the signal or not? I read in a reference that the spreaded signal appears as noise for the unintended receiver. But this only when viewed in frequency domain. Does the jammer depends on frequency domain measurements for detecting the signal, not on time domain power measurements?


texts: When negative frequencies cannot exist in practice, why do we calculate power in negative frequency?

texts: Suppose a signal with two-sided power spectral density $N/2$ is passed through an ideal LPF having cutoff frequncy  $B\,\mathrm{Hz}$. The output power is equal to $N/2 \times 2B = NB$ Watts. Why not $NB/2$?


texts: How to detect Frequency Sweep Direction when Channel is unknown

texts:
So I would like to try to modulate data using the direction of frequency sweep (a data bit of 1 would be an up sweep, and 0 a down sweep). Encoding my transmit signal is fairly easy, but i'm not sure how I would actually go about decoding the received signal. To help make things clear, please see the figure below where I have and 8 bit data stream, and the spectrogram of the transmitted and received signal. For this example my sampling frequency is $4000\textrm{ Hz}$, and my frequency sweeps are going from $0$ to $2000\textrm{ Hz}$

I've also added in a channel aspect the the figure above. My channel will be unkown, and will change from time to time (but will be constant during any single transmission). The figure below shows the channel response, which is just 2 stop bands.

In the first figure, I and can clearly see the patern of an up sweep and down sweep, but I'm not sure how I would go about actually decoding in an algorithm. I've done work with differential encoding (DPSK) and know how a method to sync and decode that kind of data, but i'm not quite sure where I would start here. The only think I could think of is to take over lapping FFT's and see which bins have the highest magnitude and try to track the direction from there.

Thanks for any input! Please let me know if I can clarify anything!

EDIT:

I tried using two matched filters as MBaz suggested such that my matched filter for my up sweep is $h$ = up sweep signal, in MATLAB it's chirp(time,startFreq,1,stopFreq). I'm then creating my matched filter output as seen below, where $y[n]$ is my output, and $x[n]$ is my received signal

$$y[n] = \sum_{i=0}^{\textrm{length}[h]} h[i]\star x[n+i]$$

I get the output as seen below, with thin spikes corresponding to the appropriate up and down sweeps. (note, it's hard to see, but there is a spike at the very beginning of the up sweep plot)

I'm hoping to get something similar to the bottom plot on the figure below, where the plot peaks at the end of the symbol, which could be used for syncing of the signal I believe.


texts: Can we choose a sampling frequency to remove unwanted noise at a specific frequency?

texts:
I am studying for my exam in signal processing. In one of the old exam papers I am told to find a sampling frequency, which will remove 80 Hz noise. The filter the exam question is based around has an impulse response $h[n]$ of
\begin{bmatrix}
-3&0&-3&0
\end{bmatrix}
My approach is to just choose the sampling frequency at a multiple of the noise frequency. In this case, I would choose it to be 320 Hz, since 4 samples gives a sum of 0. When doing a DFT with the known impulse response it would also result in a net 0 frequency response.
Does this seem like a sustainable approach?


texts: Connecting DWT matrix operations to Daubechies D4 scaling $\left(\phi(x)\right)$ and wavelet $\left(\psi(x)\right)$ functions by plotting

texts:
I'd like to duplicate the plot in the Wikipedia entry for Daubechies D4 scaling and wavelet functions:

I don't have Matlab or Mathematica, and my question is not about hacking the curve, but rather understanding the recursion that generates it. So I wonder if it'd be OK to ask for some pseudo-code implementable in Python or R, basically following the steps in the cascade algorithm, or some breakdown of the steps involved. Ideally I'd like to be able to write some ad hoc code (possibly in R, Python or Octave) to get an idea of the complexity of the process, and the difference between the differences between the discrete transform and the continuous plots above.

I see that there are initial conditions, which are resolved as an eigenvalue problem, and coefficients already worked-up for D4; yet, there is no closed formula, and an iterative numerical process is carried out.

This may involve:

Finding function $H_0,$ such that $\hat\phi(2\gamma)=H_0(\gamma)\hat \phi(\gamma),$ with $H_0\in L^2(0,1),$ and $\hat \phi(\cdot)$ standing for the Fourier transform.
Finding coefficients $c_k,$ such that

$$\begin{align}
\frac{1}{\sqrt 2} D^{-1}\phi &=\sum_{k\in\mathbb Z} c_k T_{-k}\phi\\[2ex]
\frac{1}{\sqrt 2}\frac{1}{\sqrt 2}\phi\left(\frac x 2\right) &= \sum_{k\in \mathbb Z} c_k\phi(x+k)
\end{align}$$

with the dilation operator $D_{k} f=\frac{1}{\sqrt k}f\left(\frac x k\right)$; and the translation operator $T_{k}f=f(x -k).$

Define $H_1(\gamma):= H_0\left(\gamma +\frac 1 2\right) e^{-w\pi i\gamma}$
Define the wavelet function $\psi$ by its Fourier transform

$$\hat \psi(2\gamma) = H_1(\gamma)\hat\phi(\gamma)$$

All of it rather complicated, and especially puzzling in its assumption of $\phi,$ apropos of which Strang mentions:

  We never discover the exact value $\phi(\sqrt 2)$. It is amazing to compute
  with a function we do not know — but the applications only require the
  c's.

After the comment by @Maxtron (thank you!) and using this post on SE, I was able to replicate the plots above of D4 scaling and wavelet functions - well... with some change of sign of the wavelet function:

import numpy as np
import pywt
[phi, psi, x] = pywt.Wavelet('db2').wavefun(level=12)

import pylab
pylab.plot(x, -psi,'-r', label='Wavelet fx, ' r'$\psi$')
pylab.plot(x, phi, '-b', label='Scaling fx, ' r'$\varphi$')
pylab.axhline(y=0, color='k')
pylab.legend(loc='lower right')
pylab.title('Daubechies D4')
pylab.show()


texts: If group delay is a constant for all $\omega$, Does this system have a linear phase response?

texts:
I want to know if a system has a linear phase response using frequency response of this system.

So I got phase response of it and also group delay.

and after that I got to know that if group delay is a constant for all ω, this system has a linear phase response.

so questions are

Group delay = $1/(1+\omega^2)$ . so $\omega=0$, group delay is $1$ and $\omega=1$ group delay is $1/2$. This is not a constant for all $\omega$ so I think this system has not a linear phase response. Is it correct?
when $H \left( e^{j\omega} \right)$ = $-e^{-j\omega} + 2 e^{-2j\omega} - e^{-3j\omega}$. Can I also know if this system has a linear phase response with same logic in question 1?


texts: how twiddles are calculated in the code snippet below from ti.com?

texts:
I am trying to port a C code for fft/ifft to a DSP (TMDSEVM6678LE) from ti.com. In one part of the code it tries to calculate the twiddle coefficients as you can see below; however, I can't understand why the matlab code is so different ...

I would appreciate any comments.

void tw_gen_cn (double *w, int n)
{
    int i, j, k;
    const double PI = 3.141592654;

    for (j = 1, k = 0; j <= n >> 2; j = j << 2)
    {
        for (i = 0; i < n >> 2; i += j)
        {
            w[k]     = cos (2 * PI * i / n);
            w[k + 1] = sin (2 * PI * i / n);
            w[k + 2] = cos (4 * PI * i / n);
            w[k + 3] = sin (4 * PI * i / n);
            w[k + 4] = cos (6 * PI * i / n);
            w[k + 5] = sin (6 * PI * i / n);
            k += 6;
        }
    }
}

The matlab code for twiddle coeffs. are as below:

% Program for generating n-length FFT's twiddle factor
% By: Denny Hermawanto
% Puslit Metrologi LIPI, INDONESIA
% Copyright 2015
function [real_twiddle, im_twiddle] = tw_gen_cn(fft_length)
    if nargin < 1
        [real_twiddle, im_twiddle] = tw_gen_cn(8);
        return;
    end
% fft_length = input('Enter FFT length:');
for mm = 0:1:(fft_length-1)
    theta = (-2*pi*mm*1/fft_length);
%   Twiddle factor equation
%   twiddle = exp(1i*theta);
%   Euler equation for complex exponential
%   e^(j*theta) = cos(theta) + j(sin(theta))

    twiddle(mm+1) = cos(theta) + (1i*(sin(theta)));
    real_twiddle = real(twiddle);
    real_twiddle = real_twiddle';
    im_twiddle = imag(twiddle);
    im_twiddle = im_twiddle';

end;
end


texts: Find the impulse response of an LTI system?

texts: why there is a 3/2 in the (b) problem? Does the impulse response even exist?


texts: What is the 2D version of Shift Theorem?


texts: Need to understand how wavelet is applied on Image matrix?

texts: i want to understand ,how to find out Discrete wavelet transform of image 4 by 4 size using Haar wavelet on paper ?


texts: What is the gain (advantage) of oversampling and noise shaping in D/A conversion?

texts:
It is clear that oversampling and noise shaping in A/D conversion can help in shaping
quantization noise.  But for D/A conversion, normally there is no quantization
and hence I did not understand its use.
What is the gain (advantage) of  oversampling  in D/A conversion?


texts: Detailed Analysis of Phase Delay of an Audio Signal

texts:
I write this question as a new reference point because it is difficult to find something exhaustive regarding this topic.

I'm talking about the calculus and interpretation of the phase delay of a recorded audio signal, i.e. NOT the phase response of a digital filter, as this topic is fully tackled around here.

I'm analyzing an audio signal of a hammer hit on a material, so my signal resembles a finite impulse response. I want to check if the phase can give me any physical information about my material, so I started analyzing the phase response, in terms of Phase Delay and Group Delay.

I calculate the phase delay of my audio signal as follows:

Take the whole audio file, call it $x(t)$;
calculate $X(f) = \mathfrak{F}\{x(t)\}$;
$f$ is scaled so as to correspond to Hertz scale.
calculate phase $\phi_x(f) = \arg\{(X(f)\}$, this gives me the phase shift expressed in radians;
Unwrap my phase response, $\phi_x(f) \leftarrow \phi_x(f) + n 2 \pi$. This ensures that all appropriate multiples of $2 \pi$ have been added to my phase response so that $\phi_x(f)$ is continuous (as much as possible);
calculate phase delay as: $\tau_{\phi x}(f) = - \frac{\phi_x(f)}{2 \pi f}$;

Phase delay should tell me the time delay in seconds experienced by each sinusoidal component of the input signal.

1) I know the time delay between the begin of my audio signal and the onset, and the "bias" of my phase delay (a time delay common to all frequencies) should be the same of the time delay I'm aware of. Practically, I obtain a phase delay value that is quite higher. If the time delay of the onset is about 200 samples the phase delay converted in samples approaches about 1800 samples. Why? Am I sure my scale is in seconds? It should be, because in phase delay I'm dividing radians for radians/seconds.

2) What's the meaning of negative values of my phase delay?

3) This is the plot of my phase delay:

Why higher frequencies (>8000 Hz) tend to be progressively more and more delayed? Has it a physical meaning or something concerned with the behaviour of the audio cards (non linear phase response toward high frequencies)

4) What about the frequencies components which "tend" to be not present in my audio signal (practically speaking it's impossible that in such a wideband signal a frequency amplitude can go to zero)? Is their phase calculus reliable as the most prominent frequency amplitudes?

5) If I take different windows of my signal starting after the onset, at first phase delay tends to zero, and this is logical, I'm approaching the beginning of my frequency components. When I pass with my windows the peak area of my audio signal the phase delay tends to become incomprehensible, e.g.:

I cannot understand why.

6) Last question, I need to apply windowing when I compute the frequency amplitude response of my audio signal to avoid frequency leakage and cross talk phenomena, I'm not applying it in case of the phase response computation, am i wrong with that?


texts: What is the shortest time period that can be used to find the frequency of a signal burried in noise and harmonics?

texts:
I have a continuous-time signal that is of the form
$$x(t)=A \sin(2 \pi f t+\phi_1)+B \sin(2 \pi (5f) t+\phi_2)+C \sin(2 \pi (7f) t+\phi_3)+n(t)$$
sampled at 50us where $f$ is the frequency I want to find out A,B and C are constants where $$0.15 \cdot A>B+C+|n(t)|$$. I don't need to know the value of any of the phases or amplitudes. I also know that the value of $f$ is 50Hz +/-15 Hz though i want to know it to at least +/- 1 Hz accuracy.

What ive tried so far:

A PLL-because $f$ is quite a low frequency and my frequency range is so wide and harmonics were present, the bandwidth of my loop filter was quite limited, I couldn't beat a settling time of 100ms(within 1 HZ).
DFT-Sampled over 100ms, the spectra was only really clear when $f$=50Hz.
Measuring the time between zero crossings(Though this falls apart when the amplitudes of the harmonics and noise gets quite high)

Code used for DFT

T = 0.1; % sample time 100ms
Ts = 50e-6; % Sampling period
fs = 1/Ts; % Sampling rate
t = 0:Ts:T-Ts; % Time vector of length 100ms
f1=50;
f2=35;
f3=65;
x1 = sin(f1*2*pi*t)+0.1*sin(3*f1*2*pi*t+pi/3)+0.1*sin(5*f1*2*pi*t+pi/3)+0.01*sin(50*f1*2*pi*t+pi/3);
x2 = sin(f2*2*pi*t)+0.1*sin(3*f2*2*pi*t+pi/3)+0.1*sin(5*f2*2*pi*t+pi/3)+0.01*sin(50*f2*2*pi*t+pi/3);
x3 = sin(f3*2*pi*t)+0.1*sin(3*f3*2*pi*t+pi/3)+0.1*sin(5*f3*2*pi*t+pi/3)+0.01*sin(50*f3*2*pi*t+pi/3);
N = length(t); % Number of elements in the time vector
F = fs/(N); % Frequency step for FFT display
f = (-fs/2):F:(fs/2)-F; % Frequency vector(same length as the time vector)

X=fftshift(fft(x1))/N; % FFT of the signal x1  (blue)
X2=fftshift(fft(x2))/N; % FFT of the signal x2 from task 2 (red)
X3=fftshift(fft(x3))/N; % FFT of the signal x3 from task 3 (black)

figure(2);
subplot(2,1,1), plot(t, x1, 'b', t, x2, 'r', t, x3, 'k'), axis([0 .1 -1.5 1.5]);
legend('x=50Hz','x=35Hz','x=65Hz')
subplot(2,1,2), plot(f, abs(X), 'b', f, abs(X2), 'r', f, abs(X3), 'k'), axis([-100 100 0 0.5]);


texts: Counting the number of objects intersected by the image border in MATLAB

texts:
I have an RGB image with various signs. My main goal is to count the signs that are in contact with the image borders.
Approach and problem
I started by loading the image [Fig. 1], then converted it to grayscale and applied a median filter to get rid of some noise [Fig. 2]. Then I binarized it with a threshold of 0.2, which resulted in Figure 3. At this time I got my binarized image, but the problem is that some parts that belong to the same sign are appearing in various regions, instead of only one. Now my goal is to merge the regions that belong to the same object , so then I could use bwlabel to count how many signs are in the image, and use imclearborder to get rid of the ones in the border, and use bwlabel again to get the difference between the two.
My approach was to use bwmorph, Dilate to dilate the objects and then try to fill them with imfill, holes. But the problem is that if I dilate them in a small amount[Fig. 4], the imfill doesn't seem to fill them, if I dilate them by a big amount[Fig 5] all the objects start to merge :(
Code
img=im2double(imread('image.png')); figure, imshow(img)
img_gray=rgb2gray(img); imshow(img_gray);
img_mediana=medfilt2(img_gray, [3 3]); figure, imshow(img_mediana);
img_bin=im2bw(img_mediana, 0.2); imshow(img_bin)
img_dilate=bwmorph(img_bin, 'Dilate', 10); imshow(img_dilate)
img_fill=imfill(img_dilate, 'Holes'); figure, imshow(img_fill)

Figures
Fig 1:
Fig 1 http://dl.dropbox.com/u/5272012/1.png
Fig 2:
fig 2 http://dl.dropbox.com/u/5272012/2.png
Fig 3:
fig 3 http://dl.dropbox.com/u/5272012/3.png
Fig 4:
fig 4 http://dl.dropbox.com/u/5272012/4.png
Fig 5:
fig 5 http://dl.dropbox.com/u/5272012/5.png


texts: x(t) and x(-t) are flipped (left to right) signals of each other

texts: $x(t)$ and $x(-t)$ are flipped (left to right) signals of each other.  but $x(2t-3)$ and $x(-2t+3)$ do not share that kind of relation as I thought incorrectly.  Can some one comment on this?  Is there a name for the relationship of this kind of signals?  Is it important enough to study this relation and give a name?


texts: Finding a signals complimentary sound / frequency

texts: I know about complimentary colors and how to get them, since colors are basically frequency / wavelengths.  Is it possible to find a signals complimentary, additive or subtractive frequencies in an audio file / from a signal.  I know how to use FFT to get the phase, amplitude, and frequency of a signal.  I'm just looking for a general idea of how I would go about creating similar color charts but for sound, we have a color charts of complimentary, additive or subtractive colors do we have one for the sound / audio range of humans?


texts: Why doesn't the sampling theorem work in this case?

texts: In this example, under "Discussion", we can see that the Nyquist sample rate will result in samples with the value zero only. But why is this? I thought the sampling theorem said that if you sample with the Nyquist sample rate or higher, then the original signal can be exactly recovered? But in this case it can't? So does that mean the sample theorem is not actually valid?


texts: Fourier Transform of ECG signal in Python

texts:
I have an ECG signal which I am analyzing using Python, as opposed to the mainstream MATLAB. So, I have digital form ECG in .dat file with .hea (header file). Below is the Fourier transform

The problem, as you can see, that it is not the correct Fourier transform. Below is my code.

import wfdb
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy import fftpack

cutoff = float(input("Cutoff: "))
cutoff = cutoff/(360.)
noiseAmp = float(input("Noise Amplitude: "))
programed_to = 5000

# load the aami3a ECG signal using special package wfdb for python.
record = wfdb.rdrecord('aami3a', sampto=programed_to, physical=False) # digital data, hence physical=False

# 2D numpy array storing the actual signal data. record is an object of wfdb which has many attributes.
# one of those attributes fs (sampling frequency) and d_signal which stores the amplitude at various times.
original = record.d_signal
# define time of 5000 length
t = np.arange(0, (1/record.fs)*programed_to, 1/record.fs)
# define a sinusoidal noise.
noise = noiseAmp*np.cos(2*np.pi*50*t)

# add noise to the ecg signal
record.d_signal = record.d_signal.flatten() + noise
# make the signal again a 2D array.
record.d_signal = np.reshape(record.d_signal, (-1, 1))

# define a fifth order butterworth lowpass filter with cutoff/fs/2 = 0.11 (tunable).
b, a = signal.butter(5, cutoff, 'lowpass', False)
# apply this low pass filter on the noisy ecg.
output_signal = signal.filtfilt(b, a, record.d_signal.flatten())

# plotting the FFT of the signal
# fft of the original signal.
org_fft = np.fft.fft(original)
freq = np.fft.fftfreq(original.size, 1/720.)
f1 = plt.figure()
A1 = f1.add_subplot(211)
A1.plot(freq, np.abs(org_fft))
plt.show()

I think I have taken Fourier transform correctly. Please help.


texts: Emulating guitar pickup position with delays

texts:
If I have a vibrating string, then there will be no motion at the endpoints. If I move away a bit from an endpoint then I'll get a signal, but particularly the fundamental will still be attenuated. Higher modes may suffer less or even benefit. This explains why the bridge pickup (closer to the edge) of an electric guitar suppresses the fundamental more than the neck pickup and sounds "thinner".

The amplitude of the signal obviously depends on the position of the pickup and the frequency. The relationship is some $sin(\frac{\omega x}{c})$ function.

The question is: if I have a string signal, can I emulate this position dependant behavior with delays?

Suppose the pickup is at $\frac 14$ of the string. I can then add a delay with an amplification of $-1$ which reflects (sic!) the wave bouncing back from the near edge of the string. I will get an attenuated fundamental and all that, but is it physically correct?

Why not do the same with the far end of the string and use a delay which stands for $\frac 34$ of the string. It will also attenuate the fundamental, but the notches will be at different frequencies.

Now I realize, that I am about to model the string itself. If this model was done correctly, I could excite the string and pick it up anywhere I want. But I already have a string signal (from a Karplus-Strong algorithm).

What would I have to do, to add just the effect of the pickup position to this?
What wiring of delays would produce this  $sin(\frac{\omega x}{c})$ frequency response?


texts: confused in difference between calculations of convolution and multiplication?

texts:
How are two signals multiplied? And how is it different from convolving two signals??

I have read the above question on dsp stack exchange and i also tried it in my matlab , code is as follow:

clc
clear all
close all
x=[2,4,1]
p=[5,1,8]
y1=x.*p   % result of multiplicaton
y2=conv(x,p) % result of convolution

But Matlab gives answer as given in that link but i am confused,about convolution steps, why only 5(first element of 2nd vector ) is multiplied with all there 3 element of first vector making up 3 different terms of final answer(as shown in attached photo) but 2nd element(blue under line) & 3rd element(red under line) of 2nd matrix do not repeat the same steps as followed for first element?why they are not multiplied by all elements of first vector as was the case with 5

?


texts: What is the relationship between the sigma in the Laplacian of Gaussian and the two sigmas in the Difference of Gaussians?

texts: I understand that a Laplacian-of-Gaussian filter can be approximated by a Difference-of-Gaussians filter, and that the ratio of the two sigmas for the latter should be 1:1.6 for the best approximation. However, I'm not sure how the two sigmas in the Difference of Gaussians relates to the sigma for the Laplacian of Gaussian. Is the smaller sigma in the former equal to the sigma of the latter? Is the larger sigma? Or is the relationship something else?


texts:
I asked a question earlier but I didn't get any answer for it. So now I am simplifying it: what are Cross-Spectral Density (CSD) and Power-Spectral Sensity (PSD)? What is their application? How can I get them in MATLAB?

  $$S_{kl}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_l(\omega)\} $$
  $$S_{kk}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_k(\omega)\} $$

  $S_{kl}(\omega)$ is the cross-spectral density (CSD) function between
  general signals $y_k(t)$ and $y_l(t)$, $S_{kk}(\omega)$ is the
  power-spectral density (PSD) of signal $y_k(t)$, $Y_k(\omega)$ is the
  finite Fourier transform of signal $y_k(t)$ at frequency $\omega$,
  $Y_k^*(\omega)$ is the complex conjugate of $Y_k(\omega)$, and
  $E\{\cdot\}$ is the expectation operator.

My earlier question was:
What does 'wavelet power spectrum', 'Auto-power spectrum','cross-power spectrum' means in wavelet application? I was studying about mode shape identification with wavelet method and these terms confused me.


texts: Python bandpass filter - singular matrix error

texts:
I've been trying to design a bandpass filter using scipy but I keep getting a LinAlg Singular Matrix error. I read that a singular matrix is one that is not invertable, but I'm not sure how that error is coming up and what I can do to fix it

The code takes in an EEG signal (which, in the code below, I have just replaced with an int array for testing) and filters out frequencies < 8Hz and > 12Hz (alpha band)

Can anyone shed some light on where the singular matrix error is coming from? Or alternatively, if you know of a better way to filter a signal like this I'd love to test out other options too

from scipy import signal
from scipy.signal import filter_design as fd
import matplotlib.pylab as plt

#bandpass
Wp = [8, 12]   # Cutoff frequency
Ws = [7.5, 12.5]   # Stop frequency
Rp = 1             # passband maximum loss (gpass)
As = 100              # stoppand min attenuation (gstop)

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')
w,H = signal.freqz(b,a)  # filter response
plt.plot(w,H)

t = np.linspace(1,256,256)
x = np.arange(256)
plt.plot(t,x)

y = signal.filtfilt(b,a,x)
plt.plot(t,y)


texts:
I've implemented an harmonizer, using Stefan Bernsee's code, which works great. I figured out the logic of finding the root note, and the right harmony in the scale.

The only problem is that, my upper voices sound like someone who has taken too much helium !

After some research, I found that, that was because I didn't preserve formants.
The problem is I don't know how to do that ! I've heard that I must use a Cepstrum, but that's pretty much it.

Stephan's code, modified by me is :

/* ***************** PROCESSING ******************* */
        /* this does the actual pitch shifting */
        memset(mSynMagn, 0, mParent.mFftFrameSize * sizeof(float));
        memset(mSynFreq, 0, mParent.mFftFrameSize * sizeof(float));

        for (k = 0; k <= mParent.mFftFrameSize2; k++)
        {
            index = k * pitchShift;
            if (index <= mParent.mFftFrameSize2)
            {
                mSynMagn[index] = mParent.mAnaMagn[k];
                mSynFreq[index] = mParent.mAnaFreq[k] * pitchShift;
            }
        }

I think I must mess with the SynMagn array but I don't know how.
I think calculating the Cepstrum, is just going through mAnaMag and taking the log ? (which base ?)
Then when I have that Cepstrum, what should I do with it ?
What should those line become ?

    mSynMagn[index] = mParent.mAnaMagn[k];
    mSynFreq[index] = mParent.mAnaFreq[k] * pitchShift;


texts: How to extract prosodic cues from a wav file using Python

texts:
How do I extract the following features as from wav files in Python -

Voicing rate (Voice rate) - frequency of phonemes while speaking, and represents the pace of a conversation
Number of autocorrelation peaks
Pitch, pitch bandwidth


texts: How to approximate a moving RMS without iterating over N samples for each required output?

texts:
I have a signal with a sampling rate of F Hertz.

From it I want to produce an RMS signal, also with F samples per second with a window size of N samples.

Using a standard moving RMS algorithm I must iterate over the N most recent samples F times a second resulting in F*N operations a second.

Is there a way to approximate a moving RMS signal that only uses F operations per second (or close to it)?

I know that for other processes, like iteration, there are 'leaky' algorithms that can approximate the output. Instead of a leaky iterator, is there such thing as a 'leaky RMS'?


texts:
I'm trying to match a picture saying "this is ducks, this is towers" by giving votes  from the hough transform of a picture.

I've found this article http://ttic.uchicago.edu/~smaji/papers/hough-cvpr09.pdf, but I didn't find a implemented code (and I'm trying several things).

Do you know anything similar to that, for example included on OpenCV?

Thanks.


texts:
lately i've been interested in FFT windowing, i've experimented various windows in Matlab, that made me more curious:

There are normalized parameters (CG, Scalloping loss etc...) for every window

How those parameters were measured? for what type of signals? what are the frequency, amplitude of those signal? what is the Frequency sampling used? what is the length N of the  FFT?

is there a standardized framework to use for evaluating windows.? what if i want to do an evaluation my self, what initial conditions i should choose ?


texts: Help in LBG Vector Quantization - Splitting Factor

texts:
I'm currently making a program for speech recognition. In the step of codebook generation using the LBG (Linde-Buzo-Gray) algorithm, I've read that the splitting factor $\varepsilon = 0.01$ (generally)
The splitting factor is used to split the centroid of the speech features according to the formulae
\begin{align}
Y_{n}^+ &= Y_n (1+\varepsilon)\\
Y_{n}^- &= Y_n (1-\varepsilon)
\end{align}
where $n$ is the index of the given codeword/centroid to be split and $Y_n$ is the codeword.
Also, after the codebook is generated, nearest neighbours are searched for each speech feature vector and the centroids are updated (basically clustering of features). This is done until the distortion of the codebook is less than epsilon.
Although my program seems to be working fine, I'm interested to know why the splitting factor is usually set to be 0.01.
Any help is appreciated. This is my first time working with codebooks and vector quantization.


texts: If a K-NN ($ k $ Nearest Neighbors) Algorithm Performs Very Well for Low $ k $, Can Something Be Inferred About the Data Set?

texts:
I am running kNN  on a very small data set for binary classification. Each class has 100 samples.
I am getting the best performances for $k=1$ and $k=3$. Can I deduce information about my data set from kNN performing very well for small values of $k$?


texts: Calculating the total energy of a signal

texts: I have a vector of signals and need to calculate the total energy. I have read somewhere that you can do this by summing up the entire vector and then finding the absolute value of the result. Is this correct? Could anyone suggest any tutorials, or, equations if this is the wrong method, thanks.


texts: Questions about discrete signal energy calculation

texts:
I'm looking at the definition of signal energy (e.g. Wikipedia, cnx.org). For discrete signals, it's defined like the following, where $x(n)$ holds the signal:

$
Energy = \sum_{n=-\infty}^{\infty} |x(n)|^2
$

So my questions:

For a windowed, finite signal like double signal[256], the sum is from 1 to 256 (or 0 to 255 in a program) rather than $-\infty$ to $\infty$, right? (I don't even know how I would sum over infinity.)
Why does the energy formula have the absolute-value operator $|...|$? The result of taking the absolute value is squared anyways to produce a positive value, so taking the absolute value seems to be pointless. Is it because the $x(n)$ can be complex, so the absolute value of a complex number would be the scalar from Pythagoras' theorem?


texts: Obtaining real numbers from FFT with same length as original signal

texts:
In this article the real valued time domain signal is transformed to the frequency domain to extract some features like mean or variance.
But after transform to the frequency domain I calculate average hourly values (original measurements are 1/hour), so I need the signal in the frequency domain with the same length as the signal in the time domain.
Later I need features as real numbers (e.g. mean, variance), so I also need real numbers in the frequency domain. To this end, I think they use PSD. But after using Scipy's PSD (Welch estimation) I get much shorter array.
Question:

How can I get real valued array in the frequency domain from the signal with the same length as the original signal?

How should I use FFT / real FFT / PSD to get this (what arguments' values should I use)?

Would taking absolute values of the FFT output make sense?


texts: Implementing rotation in frequency domain and map it back to spatial domain

texts:
Please consider the following small example:

x=im2double(imread('cameraman.tif'));
X=fft2(x);

Xr=imrotate(fftshift(X),90);

xRec=ifft2(ifftshift(Xr));

figure;
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
imshow(xRec);
subplot(2,2,3);
imshow(log(fftshift(abs(X))+1),[]);
subplot(2,2,4);
imshow(fftshift(log(abs(fft2(xRec))+1)),[]);

%magnitudes of both have same total energy
sum(sum(abs(X)))-sum(sum(abs(Xr)))

the result xRec is indeed rotated by 90 degrees when reconstructed from rotated Fourier domain, but a dark band appears in the center. Can someone help me to implemented in the correct manner? The goal is of course to show the Fourier rotation property. That rotation in one domain is equivalent to rotation in the other.

even more interesting, if i remove the ifftshift I almost get my rotated cameraman but with artefacts !

also some trick with polar coordinates is mentioned here (slide 22) ... but I would need help ...

https://www.kth.se/social/upload/528393f5f276545b305c8dfd/13_lecture05b.pdf


texts: Resources for learning about filters from scratch?

texts:
Topic of filters is very important in signal processing
So i want to learn about filters from scratch and then eventually reaching FIR and IIR filters
I have some idea about low pass,high pass,band pass,band stop types but often get confused how they relate to butterworth, chebyshev and other such types?
Apart from wikipedia, please recommend books chapters or websites for learning about filters


texts: I have done with iris localization. Now I need to perform iris normalization, but I don't have any idea about it. How to implement iris normalization please?


texts: How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?

texts:
I would like to calculate a delay between two signals in MATLAB. Normally, I would do this with a cross-correlation (xcorr), and then calculate the position of the correlation peak, but in this case I know that the delay is smaller than one sampling period. Therefore, the method I propose will yield either $0$ or $1$ sample as a delay.

One solution I can come up with, is to construct an interpolated version of the two signals, and then calculate the delay. However, then I still rely on the sampling period of this upsampled signals. Another possible solution is to calculate the cross-correlation peak for different blocks, and then calculate the mean delay. However, I am not sure whether this solution will converge to the right delay. A third solution I can come up with, is to apply a fractional delay to one of the two signals by means of a sinc reconstruction, and then calculate for which fractional delay the correlation peak is the largest. This requires a lot of calculations however, and I still rely on the fractional delays I apply (so I end up with the same problem as with the upsampling).

Therefore, I am wondering whether there is a more elegant or precise solution?


texts: Obtaining the size of objects from an image

texts:
My objective is to be able to obtain the size, in pixels, of objects in images I'm receiving from a camera.

I've tried two approaches thus far, both working to some extent, and I'm wondering if there's any other approach you can suggest.

The first approach was using the CIE Delta E formulas. I would parse through the BGR 24 image data and then locate the appropriate pixels that form part of the object I wish to size. This approach works well when I know the object colors and positions beforehand, not so much when I don't. Another drawback of this approach is that calculating the Delta E for every pixel is extremely computationally expensive.

To circumvent my reliance on having to know the colors and positions beforehand, I used Connected-component Labeling to find the objects and size them.

Here's an example of what I'm doing using connected component labeling.

As you can see from the image above, I'm making four measurements, horizontal, vertical, left diagonal, and right diagonal. After I've sized the object I know where it is in the image, I then debayer the image into a BGR24 image and take several color samples. Debayering the image constitutes about 50% of the processing time, which is more than I prefer.

Now, my question is: What other approaches are there for finding the size in pixels of objects in an image?


texts: How to scale measured QPSK points to determine magnitude error?

texts:
I am trying to determine the magnitude error between demodulated I,Q values and reference values on a constellation plot for a QPSK signal. The reference points are at (1,0),(0,1),(0,-1), and (-1,0). Of course as the amplitude of the input signal increases, the measured points move further away from the references. My (limited) understanding is that measured points would be near the expected to make the resulting constellation useful - where the proximity is affected by the quality of the signal, but not the amplitude.

Is there a standard scaling method used to shift the measured output to the range of the references for this purpose?


texts:
I'm confused with the last step of extraction MFCC where we have to take DCT or IDCT.

Cepstrum is defined as IDFT of log(power spectrum). Since MFCC are cepstral coefficients so I think we have to take IDCT at the last step.

It's confirmed by several articles. For example this one (point 2.4 we see F^(-1)) and this one (scheme on the second page). In fact in first link formula below F(^-1) is DCT-2. Why? As wiki says DCT-3 is an inverse transform.

And if we look at implementations they take DCT-2 too. I checked two of them: pyAudioAnalysis (stMFCC method) and implementation (mfcc method) by the guy who writes for practicalcryptography.com.
And this is also confirmed by several links: here and here (second page).

Well.. who is right?


texts:
I have a Biquad II filter implementation in my embedded environment and I'm trying to find some equations to generate the coefficients for it.

I have already found the Audio Eq Cookbook, which simply describes how to calculate all of the coefficients for Type I.

LPF: (Taken from Audio Eq Cookbook, Type I Equations)

omega = 2*PI*frequency/sampleRate

sin   = sin(omega)
cos   = cos(omega)

alpha = sin/(2*Q)                                     (if Q is specified)
      = sin*sinh[ ln(2)/2 * bandwidth * omega/sin ]   (if bandwidth is specified)

            b0 =  (1 - cos)/2
            b1 =   1 - cos
            b2 =  (1 - cos)/2
            a0 =   1 + alpha
            a1 =  -2*cos
            a2 =   1 - alpha

What are the equations for type II coefficients?